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ABSTRACT 


The propagation of the charge-exchange plasma for an electrostatic 
ion thruster is crucial in determining the interaction of that plasma 
with the associated spacecraft. A model that describes this plasma and 
its propagation is described, together with a computer code based on 
this model. 

The structure and calling sequence of the code, named PLASIM, is 
described. An explanation of the program's input and output is included, 
together with samples of both. The code is written in ANSI Standard 
Fortran IV. 
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INTRODUCTION 


Ion thrusters can bo used in a variety of primary and ouxiliary 
spoce-propulsion applications. A thruster produces a charge-exchange 
plasma which can inter&ct with various systems of the spacecraft. In 
order to understand those possible interactions, o detailed knowledge 
of the plasma propagation la required. 

The production of charge-exchange ions by thrusters has been 

1 

understood for some time. Fast ions from the thruster interact with 
slow neutrals that are also escaping, resulting in the production of 
ions that initially have only a thermal velocity. The electric fields 
within the ion beam cause these ions to move approximately radially out 
of the ion beam. These charge-exchange ions leave the ion beam along 
with electrons supplied by the neutralizer, the combination constituting 
the charge-exchange plasma. The propagation of the charge-exchange 
plasma depends on several factors, Including the initial thermal energy 
of the ions, the distribution of ion production along the beam, and the 
potentials and geometry of neighboring spacecraft surfaces. 

In the THEORY section of this report, the geometry of an idealized 
spacecraft with an ion thruster is described, together with the simpli- 
fications and definitions used in modeling the ion beam. The distribution 
function used for charge-exchange ion production is also presented, 
along with the barometric equation that relates the variation in plasma 
potential to the variation in plasma density. The numerical methods and 
approximations used for the calculations are then discussed. This 
section describes the main calculation subroutine, CALC, and the dis- 
placement calculatl-on subroutine, CALCD. 
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In fcho PROGRAM STRUCTURE section, a flowchart is provided that 
diagrams the calling soqucnco of the modules; also presented are 
detailed descriptions of each of the modules. A guide to using the 
program is presented in the INPUT AND OUTPUT section of ihia report. 
Descriptions of the calculated results are presented in the INPUT AND 
OUTPUT and VERIFICATION sections. 

Also presented is a method of obtaining better resolution in the 

upstream region. The high~resolution option of the program simulates 

only the upstream region. This option utilises previously calculated 

« 

trajectories as boundories for the region to be simulated at higher 
resolution (see notes in the computer code). 

The VERIFICATION section of this report also compares experimental 
and analytic results with those obtained by the corapuLer code. 

Factors limiting simulation accuracy are also discussed. 

An analytic solution is derived for the case of an infinitely 
long cylindrical beam with a uniform distribution of charge-exchange ion 
production along the beam. Expressions are obtained for the radial 
variations in ion density and velocity, permitting a direct comparison 
with results from the computer code. This analytic solution is 
described in APPENDIX A and used in the comparison described above. 

It should be noted that this final report provides a complete 

2 3 

description of the program and supersedes previous reports, * All of 
the Information necessary to use the program is contained herein. A 
glossary of the variables used in the computer code is provided in 
APPENDIX C and the computer code is listed in APPENDIX D, 


ORIGINAL* PAGE IS 
dF POOR QUALITY 



3 


THEORY 

The Intciraetion of an Ion thruscet with other eoroponcnto of an 
elecU’icoXXy propoXXed opacoeroft through the pXnoma ourrounding a 
opacccraft has been otuv'.ted for oome tlmoi The tranoport of eXectrono 
from the ion beam to a soXar*“array ourfacc wao treated firot by Knauer, 
ct ax/* ao on cXcctron spaeo-eharge-fXow problem. Mcaoured cXcctron 
eurrento, though, were found to be much higher than eaXeulated by 
Knauer. The difference wao due to the presence of a charge-exchange 
pXasma. 

Charge-exchange ions are produced when fast beam ions pass near 
oXow escaping neutraXs. The fast neutraXs that result usually present 
no problem, and escape following the directions they had as ions. The 
slow charge-exchange ions that are produced, though, initially have only 
the velocity of the thermal neutralo. Small electric fields within the 
ion beam result in the charge-exchange ions leaving the beam in approxi- 
mately radial directions. These charge-exchange ions, together with 
some escaping eleetrono, form the charge-exchange plasma that surrounds 
an electrically propelled spacecraft, 

The production rate for the charge-exchange ions was first calcu- 
1 

lated by Staggs, et al. The capability of the charge-exchange plasma 

to transport electrons to other parts of the spacecraft was experimentally 

5 

evaluated by Worlock. at al. Some detailed trajectories of charge- 

A 

exchange ions have been examined by Komatsu, et al, Experimental studios 

of the charge-exchange plasma distribution, particularly upstream of the 

ion-beam direction, have been conducted by Kaufman, and Carruth, et 
10-12 

al. Several studies included a correlation of plasma properties in 

terms of the distance from the thruster and the angle relative to the 
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beam diraation. Theoretical ctudioo o£ th& charge exchange plasma 
have been carried out by Rubinoon, ct and Katz, et al,^^ The 

latter treat the lono, numerically, ae o cold fluid in contrast to the 
use of calculated ion trajectories and density gradients* 

The physieal proceoaeo Involved in the ehargc-exchange plasma have 
become well understood ao a result of the various studios that have 
been conducted, The electron population outside of the beam agrees with 
the "barometric" equation 

"o “ %,rrf “P(-lV/k*e) 0) 

15 

which wan introduced by Sellen, ct al. and verified by Ogawa, ct 
10,4./ population within the beam and by Kaufman for the 

population in the charge-exchange plasma. The plasma potential V in Eq. 
(1) is taken to be zero at the reference electron density n^ The 

electron temperature T^ in the charge-exchange plasma has been found to 

7 

be about half the value in the ion beam. The electron temperature in the 
ion beam varies with thruster size and ranges from about 7 eV for a 5 cm 
thruster to 5 eV for IS cm and 0.35 eV for 30 cm. Also, q is the elem- 
entary unit of charge and k is Boltzman’s constant. 

The experimental validity of Eq. (1) is consistent with the low 
density and long mean-free paths in the charge-exchange plasma. The 
decreasing plasma density with increasing distance from the thruster 
forms a potential well for the electrons, so tlut many transits of this 
region are probable before an electron escapes. The many transits 
permit randomization of the electron population to a single Maxwellian 
distribution by Coulomb collisions. 
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The ol the ehatge-oxehange pAastna io Jorge compared to th?i 

Debye ohieAding diatonee, which means that the electron density must 
everywhere be equal to the ion density • Inasmuch as the Ions only move 
outward £rom the thruster, ttelr motion is cosentlally collioionleos and 
governed by the potential distribution £rom Eq. (1). 

The approach used in this study has been to assume a cylindrical, 
axially symmetric ion beam, with the charge-exchange ions leaving the 
beam with a uniform velocity in the radial direction. The coordinates 
and simulation boundaries are defined in Elg, 1. The current density of 
these chargo-exchange ions at the cylindrical beam boundary io a function 
of the distance downstream from the thruster. The total ehargo-exchange 
current is distributed among the total number of trajectories, with 
this total number specified as an input parameter, N. Approximately 
fifty percent of the charge-exchange ions are generated within one beam 
radius of the downstream end of the thruster, so about half of the 
specified trajectories will initially start in this region. 

A trajectory represents the path of a single representative ion on 
which acceleration is produced by electric fields in the plasma. These 
fields correspond to potential gradients induced by gradients in the 
plasma density, as indicated by Eq. (1). Density gradients are used in 
two separate calculations. The component of the density gradient along 
the path provides a potential gradient which serves to change the ion 
velocity in that direction, while the component of the density gradient 
normal to the path provides a potential gradient which modifies the 
direction of the path. The forces on the ions acting normal and parallel 
to the path are resolved into x and z components to obtain the resultant 
force acting m the ion. 
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In the simulation used heroin, the ion path is represented with a 
stepwise progression away from the beam and the trajectories are advanced 
from left to right starting from the end of the thruster. It is assumed 
that, with small enough step siaies, following the Ions through one pass 
of calculations is sufficient. (The validity of this assumption is 
partially checked later by comparison with eKperimental and analytic 
results.) From a physical viewpoint, the ions are moving at, or above, 
acoustic velocity, so disturbances should not propagate in the upstream 
direction. Also, the extent of the plasma is very large compared to the 
Debye shielding distance, he>ce electric fields at the flow boundaries 
should not extend into the bulk of the plasma. The distances to neigh- 
boring paths are Important parameters, in that they are used to determine 
densities. As indicated In Fig. 2, a normal to the path being incremented 
Is extended in both directions. This nermal is used to calculate the 
distances to neighboring paths, on the right and left of the path being 
incremented. In this report, right and left are defined in terms of 
relationships upon leaving the ion beam, with the viewing direction in 
the direction of charge-exchange ion motion. The distances to the 
neighboring paths are obtained by calculating the distance from the path 
being Incremented, along the normal to this path, to the neighboring 
path where the normal intersects it. If the neighboring path was not 
Intersected by a normal, as in the right side of Pig. 2a, then the 
neighboring path is extended linearly from the last Interval and the 
distance is calculated. If the normal intersects a neighboring path 
below the final two ion positions on that path, as is the case for the 
left hand paths in Figs. 2a and 2b, back stepping is used. This allows 
the distance to be calculated to the line the normal actually intersects. 
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(a) Linear extrapsOlation used on right, bock stepping 
used on ieft. 



(b) Paths intersected both sides { path iterated again ), 
back stepping used on left. 


Fig. 2. Geometry for evaluations of distances between paths. 
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If Interoections arc found on both the right and left, as In Fig. 2b 
(linear cKtrapolation not used) , the path currently being iterated will 
bo iterated again so it can keep in stop with its neighbors. Finally, 
if the normal intersects the neighboring path between the final two ion 
positions, as in the right hand path of Fig. 2b, the intersection point 
is used as it is. 

The density is inversely proportional to the distance between neigh- 
boring paths, the ion velocity and the radial distance. The latter rela- 
tionship is duo to the axial symmetry and the use of only one trajectory 
for each axial location. The density or the left is thus given by 

“ C/Adj^xv^, (2) 

where C is a constant depending on operating conditions and the number 
of trajectories specified, Ad ^s the distance betvteen the path being 
incremented and the path on the left, x is the radius and v^ is the ion 
velocity. The density on the right is defined in a similar manner, except 
that and Ad^^ are used. After the radial distance, x, and the ion 
velocity, v^, are calculated, the densities to the right and left are cal- 
culated using Eq, (2). The two quantities, n. and n„, are averaged to 
get the plasma density at the point under consideration. 

The plasma potentials to the right and left are then calculated 
using Eq. (1). Here n^ represents the density just calculated to tlie 
right or left and n^ represents the initial density to the right 
or loft of the path being incremented. The force normal to the path 
direction is then, 


JL 


-qAVj^/Adg 


(3) 
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wUeve io the difference between the potentials on the right and left 
and u,d , is the smaller of .^(t and Ad^. This choice for Ad„ was included 
to accommodate radical changes in the displacements arising when the 
perpendicular displacement intersects a boundary. The effect of this 
choice as compared to averaging the displacements was found to have a 
negligible effect on interior paths. This force is resolved into x avid 
2 components. 

The force acting parallel to the path can be calculated in a fashion 
similar to Eq. (3), 

F|j - -q^V||/Adp (4) 

vdvere ^sV y is obtained through Eq, (1) with potentials being those at the 
point presently under consideration and the previous point and the 
distance between these two points, This force is also resolved into 
X and z components. 

The normal and parallel forces can also be set equal to the rate of 
change of momentum: 

"i,r 

with m the ion mass, 4vj^ U the velocity component generated normal 
or parallel to the path direction, and 4t the size of the time interval 
list'd In the iteration. Equating these two force expressions yields 

These are the velocity changes normal and parallel to the path direction 
for the present iceration. Similar expressions hold for the x and z 
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components of the velocity, and Av^, where „ ia used instead of 
Fj^ II* total velocity and its components are calculated from Av^^ and 

Av^. The new ion position is tljsn calculated, using linear expressions, 
thereby incrementing the path* 

It was necessary to consider several special cases in the execution 
of the displacement algorithm. Three were mentioned in reference to 
Figs, 2a and 2b in discussing the calculation of distances to neighboring 
paths. Those were linear extrapolation, back stepping and iterating a 
path again if intersections were found on both sides. Other cases 
involve the extreme right and left trajectories that have not inter- 
cepted a boundary, Without special treatment, these cases would result 
in an undefined density on one side of the path because the normal will 
intercept a boundary instead of another path. The boundary is treated 
as another path with one exception. If at any time a path would be 

repelled by the boundary, the direction is left unchanged. This approx- 

* 

imates a plasma sheath which would be present at such a boundary. 

In general, both the distance between trajectories and the distance 
from a trajectory to a boundary will be much larger than the Debye 
distance. The accuracy of the simulation should therefore be considered 
questionable at any location where the distance between trajectories 
approaches the distance to a boundary. A better approximation in such a 
location might be obtained by extrapolating from deeper within the 
charge-exchange plasma. It would also be possible to use more trajec- 
tories, so that the space between them would be reduced. 

The distribution of charge-exchange ion production along the axis 
ia assumed to be proportional to the neutral density on the axis. This 

neutral density for a single thruster (no overlap of neutral effluxes 

7 8 

from adjacent thrusters) is ’ (see Appendix B) 
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n(z) 


( 3 ^ 4 . r , 2)^/2 


( 7 ) 


whoro is the beam radius and is a constant for a given combination 
of beam diameter and neutral loss rate. This function decreases rapidly 
with increasing z, due to the rapid divergence of neutral atom paths in 
free molecular flow. The beam radius, is an important parameter in 
this simulation, because approximately half of the total charge-exchange 
production occurs within about one beam radius of the thruster. This 
means that half of the charge-exchange trajectories will originate 
within the some distance. 

In determining the locations for the origin of ion trajectories 
along the axis, the integral of Eq. (7) is used 


00 

/ "^2)dz - n^r^ . 
o 


( 8 ) 


The region simulated is finite, so that not all of the Integral can be 
represented. The region to be simulated was defined so that 95 percent 
of n^r^ is contained within this region. For N trajectories making up 
0.95 with each trajectory located at the median of the density 

interval that it represents, the following expression holds 

2N ^i+1 ^1+1 

■| S / n(z)dz « N r n(z)dz » 0,95 n r. /2 . (9) 

i-O 3^ 

For the first trajectory, the starting point is at the right end of the 
thruster (left end of the ion beam, see Fig, 1). To calculate each 
successive z, the expression used is 


0.95r|^/2N = - 3^ + • (10) 
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For the first trajectory, i » 0 and i + 1 1, The value of is the 

right end of the thruster and the first trajectory is started at 
The second trajectory is started at z<^, third at and so forth. 

The initial velocity upon leaving the ion beam Is the Bohm velocity, 

Vj - . (11) 


where T is the electron temperature in the ion beam and m. is the ion 
mass. The constant, C, in Eq, (2) is obtained using the following 
procedure. The total production rate of charge-exchange ions, for a 
uniform beam current density profile, is given by” 



2 - 


« 2J, (1-n )a /irr. n,.q v 

b u ce' b i 


( 12 ) 


where is the ion-beam current (A) , is the propellant utilization, 
0 ^^ is the charge-exchange cross section (m ) , r^ is the beam radius (m) , 
q is the absolute electronic charge (C), and v^ is the mean neutral 
thermal velocity, (8kT^/Trm^) ' (ra/sec) . With typical values for Hg 
neutrals and singly charged ions used, 

= 6.18 X 10^^ 

The charge-exchange cross section usually decreases with ion energy. 

+ 

The value used for Eq. (13) corresponds to Hg ions at about 1,000 eV. 

The plasma density can be related to this production rate by 

n = N /27rAd xv.N , (14) 

ce ml 
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where is the local mean spacing between trajectories (m) , x is the 
radius (m), is the local ion velocity (m/scc) and Vi is the number of 
trajectories aimulatcd. Substituting for 

\ _i_. 

where the quantity enclosed in the paronthosos is the constant C» 
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PROGRAM STRUCTURE 

The simulation is performed by a eoroputor program written in 
otondard Fortran IV that is listed in APPENDIX D. The main driver 
program, PLASIM, establishes the oequenelng of computations in the 
simulation and makes calls to the major subroutines vfhieh perform the 
required naleulations. This overall flow is illustrated with a flow 
ohart in Fig. 3. 

First, all of the fixed parameters are stored in the COMMON bloeks 
by the subprogram BLOCK DATA. The control parameters which define the 
type of computation, disposition of results and termination, along with 
physical data used in the simulation, such as the thruster and accel- 
erator system dimensions and the plasma characteristics, are input 
through subroutine READER. The parameters and data are read from card 
Images. 

For a typical simulation, the next call is f;o subroutine INIT, 
which initializes the constants and arrays, calculates the coordinates 
of ion trajectories (paths) at the beam edge and performs the first 
iteration, thus calculating the second set of coordinates on the paths, 
Calls are then made to subroutine WRIT to output the heading, thruster 
schematic, initial parameters and results of the initialization and 
first iteration. 

PLASIM next begins the staging and successive iterations by calls 
to subroutine CALC. Subroutine WRIT Is also called to output a heading 
for information on the progress of the computation. After the completion 
of a specified number of iterations, termed a stage, the contents of the 
coordinate and density arrays are written to an external mass storage 
file, PATHS. The arrays are reinitialized and the next stage commenced. 
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Fig. 3. Calling sequence for program PLASIM. 
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After completion of oil the etoges, the plotting routinee are called 
according to the status of the control parameters set at the beginning# 
Another call to READER is made to determine the next actipni either 
another aimulation run or a normal termination. 

The main calculation aubroutine, CALC, sets up an accounting 
procedure for the paths, the iteration number, the activity of the 
path, and the boundary condition, in addition to carrying out calcula- 
tions leading to new ion positions. Displacements from the current 
path to adjacent paths are returned by a call to subroutine CALCD and 
then the densities on both sides of the present path are calculated. 

The potentials, on both sides of the path and at the present and prior 
ion locations, are obtained, followed by the calculation of forces 
acting perpendicular and parallel to the path, respectively. The forces 
are resolved into z and x components to calculate the new z and x 
components of the ion velocities. A new ion location is computed from 
the z and x velocities and a call is made to subroutine BOUND to ascertain 
if it is within the simulation boundaries. 

Two plotting methods are included, one for a line printer, LNPLT, 

and one for a Versatec electrostatic plotter, VRSPL. The subroutines 

18 

called in VRSPL are described in the appropriate literature. These 
may be not only device-dependent, but site-specific as well, They are 
included to indicate a possible preparation procedure to use in presenting 
the- simulation results in graphical form, 

There are several error detection segments in the code which output 
informative messages whenever unrealistic values appear in certain 
variables (see subroutine BLOCK DATA) . 
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Thuro ai'o three typos of simulations defined by Che control parameter 
KEY* KF<y*"“l nonoratos a normal simulation with s non --uniform initial 
dennity diotribufelon, KEY*»“1 Konerates a elmulaelon with a uniform 
initial density distribution and KBY^O generates an upstream simulation 
which utillacD a given path from a previous run as a boundary and 
simulates paths upstream from that path. A call with KEY*0 causes a 
normal termination including normal clearing of the plotting buffers, 

The following oubseetions discuss the various subroutines in detail# 

PLASIM 

This is the main driver routine for the simulation of a charge- 
exchange plasma surrounding a brood beam ion thruster. This routine 
controls and sequences the computation and lefines the program structure 
and staging by calling various subroutines. It is outlined in Fig. 4, 

The first statement (two lines) of the program is a non-ANSI statement 
which declares the files required for input and output. This will bo 
converted to a comment statement to prevent errors because of non- 
standard FORTRAN. 

READER: The first call Is to subroutine READER which reads in the 

run specifications and Information, INIT: The next call is to sub- 

routine INIT which initializes various parameters and performs the first 
iteration. CALC: The next call is to subroutine CALC which computes 

the ion positions along the trajectories. CALC is called NUMIT times. 

This completes the first stage, (NUMIT is the maximum number of itera- 
tions to be performed on any given path during any particular stage,) 
Information for the first (NUMIT - 10) iterations is written to an 
external file, PATHS. Gore memory is then reinitialized by transforing 
the results of the final ten iterations to the core space for the 



Fig. A. Flow chart of main driver routine, PLASIM, 
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Fig. k. Continued* PI.AS1M. 
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initial ten. The results of these ten iterations aet os a base for 
another set of (NUMIT - 10) iterations. CALC is called another (NUMIX - 
10) times to fill core storage with another set of PATH coordinates 
which completes the second stage. Then another set of (NUMIT - 10) 
results are written to PATHSi core is reinitialized again and another 
stage is run. This is done until all paths intersect boundaries, 
Finally, the plot indicator, ICLPLT, is checked to determine which 
plotting routine to use to output a plot of the path coordinates, if 
any, 

Two blank cards are always placed at the end of the data deck to 
signal a stop, The word size dependent variables are IFl, IF2, INFO, 
ITITL, and IW, See BLOCK DATA for instructions on word size dependent 
variables and code generated error messages, 

BLOCK DATA 

This subprogram loads data into labeled common storage at compile 
time through the data statements. Error message information and the 
Instructions on word size dependence have been Included. This depen- 
dence was motivated by the requirement for transfer of alphanumeric data 
for labels to the plotting subroutines of Versatec and IMSL. The basic 
unit for input was chosen to be 80 characters, a full data card. These 
80 characters must be packed cotitint?oui?ly into words of an array which 
will be transferred to the exiT^rual plotting subroutines. The arrange- 
ment described below accomplishes this but requires a change in the DATA 
statements for computers with word size not equal to 10 characters. 

The word size dependent variables are; IFl, IF2, IW, INFO, and 
ITITL, IW is the number of words required to generate ^0 characters. 

To convert to a machine using different word size, modify only IW, IFl 
and IF2 in first two data statements below. The third data statement 
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io £oi‘ I/O l)u£f(5i’H and the fourth data statement is for values of 
comitanto. 

DATA IW, mO), m(2) /8, 6HC8A10), 111 / 

DATAIF2(1), IF2(2), m(3), IF2(4) /10ll(8A10/2(AA,4H10)) ,1H ,1H / 
DATA IN, lODT, IFATHS /5, 6, 7/ 

DATA BK, Q, PI /1.3806E-23, 1.602E-19, 3.141593/ 

When a code generated error moaaage la called, the first line of 
the error output will bo of the form, 

ERROR NNN 

where ’NNN* ia one of the following integers, 

207 See Subroutine WRIT 

410 ■''ee Function Subroutine DS 

412 See f'unction Subroutine DS 

521 See Subroutine CALGD 

522 See Subroutine CAbCD 

523 Bee Subroutine CALGD 

524 See Subroutine CALGD 

525 See Subroutine GALCD 

526 See Subroutine CALGD 

527 See Subroutine CALC 

528 See Subroutine CALGD 

529 See Subroutine CALGD 

530 See Subroutine CALC 

610 Sec Subroutine BOUND 

612 See Subroutine BOUND 

and the referenced aubroutine is the colling routine, 


READER 

See INPUT AND OUTPUT aecUion. 


INIT 

This subroutine initializes the necessary variables and performs 
the first iteration, It is outlined in Pig, 5. The constants, which 
include the Bohm velocity, velocity of the thermal neutrals and the 
stop sire, are defined first. Then the z and x coordinates ot the ion 
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Fig, 5. Flow chart of initialization routine, INIT. 
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Fig. 5. Continued, INIT. 
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trajectory exit points are obtained and ZBOUND is defined. If this is a 
high resolution upstream run (KEY>0), ZBOUND is redefined. If this is a 
uniform density run (KEY*-1) , the ion trajectory exit points are redefined 
to be unif'^rmly spaced. Following this, the velocity arrays and iteration 
counter are initialized, the second ion positions on each path are calculated 
and all the paths are set to active. The initial densities between the 
paths and the initial densities on the paths are then calculated. Finally, 
the heading, thruster schematic, initial parameters and results of the init- 
ialization and first iteration are output, then control is returned to PLASIM 

CALC 

This is the main calculation routine. It is outlined in Fig. 6. 

This subroutine uses the arrays in blank common along with subroutine 
GALGD to determine the next position of the ion being considered. The 
densities and potentials to the right and left of the present path are 
calculated first. Then the force acting perpendicular to the present 
path is calculated. If the boundaries repel the path, the perpendicular 
force is set equal to zero. The potentials and forces acting parallel 
to the present path are obtained next. The total force (sum of perpen- 
dicular and parallel components) acting on the ion is calculated and then 
the velocity components along the x and z axes are obtained. The speed on 
path one is normalized to 1, 2 times the speed on path two when the speed 
on path one is 20% greater than the speed on path two. Finally the next 
ion position is calculated. Subroutine BOUND is called to make sure the 
new ion position is Inside the boundaries. The results are printed every 
IGLWRT iterations. IFLAGl is checked (see GALGD) to see if intersections 
were found to both the right and the left. If intersections were found 
to both sides, the path is iterated again. 





Fig. 6. Flow chart of main calculation routine, CALC. 
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Fig* 6. Continued, CALC* 
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CALCD 

Thin subroutine, which is outlined in Fig, 7, computes the perpen- 
dicular displacement from the present path to the adjacent paths » The 
flagn arc defined in this routine. 

First, the flags are Initialized and then the slope of the line 
perpendicular to the present path, at its end point, is calculated, If 
the path under consideration is the right- or left-most active path, 
function subroutine DS is used to get the displacement to the boundary. 

When on interior path is under consideration, displacements are 
obtained to the right- and left-hand paths. Here, the slope of the 
adjacent path, between the last two ion positions, is calculated. The 
intersection point between the line perpendicular to the path under 
consideration and the line formed by the last two points of the adjacent 
pa.th, is calculated. Tests are run to determine if linear extrapolation 
(adjacent path extended to find an intersection point) or back stepping 
(Interpolation between points on the adjacent path previous to the last 
two points) is to be used. The flags are set, the displacement of the 
adjacent path is calculated and then control is returned to subroutine 
CALC. This subroutine returns the values of the flags, the displacements 
to the right and left and the angle that the perpendicular makes with 
the horizontal, to subroutine CALC. 

BOUND 

This subroutine checks the point (z,x) to see if it lies inside the 
defined boundaries of the simulation. If (z,x) does lie inside the 
defined boundaries, no changes are made and control is returned to sub- 
routine CALC. If the point lies outside the boundaries, the path status 
is set equal to the iteration number. In addition, if (z,x) lies on the 
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. 7. Flow chart of displacement calculation routine, CALCD. 
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Fig, 7. Continued, CAtCD. 
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Fig. 7. Continued, CALCD. 
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first or laofe active path ond lies outside the boundaries, no changes 
are made ond control io returned to subroutine CAhC. 


WRIT 

Thio io the main output routine. It consists of five subsections 
with the value of the passed parameter, KE, determining which section is 
called. KE«1{ the first three pages of output arc generated, consisting 
of the beading, a thruster schematic ond the constants and input data, 
KE«2; the interim status of the major variables is output from subroutine 
INlT. KEa3s the heading is printed for the results of one pass of cal- 
culations. KE«^: a file of path coordinateo is created to be used in a 
high resolution upstream run, KE°5! a file of position-density triplets 

^n »*n ueorl fnv tana TMniW AMn rtti'T'DJl'i’ 


DS 

This function subroutine finds the perpendicular displacement from 
the first or last active path to the boundary, CALCD constructs a 
perpendicular to the present path at the present point (z,k) with a 
slope of SLOPEP. Function DS then extrapolates this perpendicular of 
slope SLOPEP to the left or right (depending on whether the first or 
last a'tive path is considered) and finds the z and x intercepts (ZINT 
and XINT) along a boundary line, ZINT and XINT are checked to see if 
they lie on or between the boundary endpoints on the boundary line. If 
they do, the displacement is calculated, if they do not, ZINT and XINT 
are calculated along the next boundary line and tasted again, This 
continues until adequate intersection points are found or all boundaries 
have been considered in which case an error message is output, The 
perpendicular displacement is returned to subroutine CALCD. 
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VRSPL and VRSPLT 

Tlicot^ oubroutineo net up the data in XION and XION for tranofcr to 
the utility plotting paekago, Veroatec, They aloo draw a oehematlc of 
the thruQtor and beami habolo for the graph and ito axes are aloo 
tranof erred and involve the word-oize dependent varlabloo. VRSPL plo4o 
from the external file, PATHS, whorcao VRSPLT ploto from eore memory. 
Thcoo routinoo are device dependent and may require considerable modifi- 
cation at other installations. 

LNPLT and PLOTW 

LNPLT sots up the data in XION and XION for transfer to the utility 
plotting routine PLOW which utilizes the line printer. It la also 
dependent on the type of printer available and may require considerable 
modification at other inatallations. 
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INPUT AND OUTPUT 

Por each Glmulatloii run, ooven data earde (or cord images) arc read 
by READER. On the second card, if KEY»0, the reading aoqucnee is inter- 
rupted and control returns to PLASIM, whore normal termination ensuos, 
The input data cards have the following formats and parameters which are 
defined in the following subsection and in APPENDIX G. 

Cord 1, PORIIAT (see the first data card in BLOCK DATA) 

Content; Description of the run, up to 80 characters. 

Card 2, FORMAT (7110) 

Content; NUMION, NUMIT, KEY, ICLWRT, IGLPLT, NTOTST, IGLERR. 

Card 3, FORMAT (6F10.5) 

Content; RB, RBOUND, RT, THRLEN, BMGLR, UTIL. 

Card 4, FORMAT (5E10.3) 

Content; TELIN, TELOUT, TTHNEU, GEXSEC, UMSION. 

Card 5, FORMAT (3F10.3) 

Content; TIMEMU, XVELMU, EVELMU. 

Card 6, FORMAT (see the second data card in BLOCK DATA) 

Content; Label for graphical output. 

Card 7, FORMAT (see the second data card in BLOCK DAT.A) 

Content; Axis labels for graphical output. 

Due to the dimensions in the COMMON BLOCK, the values of the 
variables NUMION and NUMIT should not exceed 41 and 150, respectively* 
Also NUMIT should not be less than 11 and the total number of itera- 
tions performed should be less than 8888. The temperatures input on the 
fourth data card should be in units of eV, all other variables require 
SI units, 

The printed output begins with three pages of identification 
including a schematic of the simulation region and a statement of the 
parameters defining the conditions of the run. (See WRIT(l).) 
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The noKt two sogmente coriaiot of INTERIM STATUS » a report of the 
coordioatOG, velocity and density for each path after the initial tration 
and first iteration, os prepared by INITi (See WRIT(2).) A provision 
is made in subroutine CAtG to output tVie counters, coordinates, velocity 
and density values for each iteration so that the progress of the 
simulation may bo followed in detail. If the results of each iteration 
are not desired, the value of the input variable ICLWRT on data card 2 
should be set to a value other than one. The value chosen determines 
how many iterations must occur before those results are output. WRIT(3) 
outputs the heading for the information output in subroutine GAIG. 

Graphical output consists of an outline drawing of the thruster and 
beam boundaries and the paths generated by the computation in the form 
of dots or lines. The paths are graphs of the coordinates contained in 
the arrays ZION, XIQN. These arrays are each a doviblo-subscripted array 
in which the first subscript identifies the ion path number, and the 
second one identifies the iteration. The density array is also available 
for plotting, but provision for such has not been Installed. 

An appropriate relative scale in terms of the number of paths and 
the iteration step size must be established to properly model the upstream 
regions. Very small steps compared to the path spacing are inappropriate 
as are very large steps compared to the path spacing. The scaling ratio 
S “ ad/Vj^^At compares the path spacing to the axial step size<> Appropriate 
values of this ratio, used herein, for the various thruster sizes arej 
5 cm thrusters S =» 2.9, 15 cm thruster; S 3.9, 30 cm thruster; S “ 3.6. 
These values were obtained by taking Ad as the separation between the 
14th and 15th paths, along the beam edge, Vg as the Bohm velocity and At 
as the time step for each configuration given in the SAMPLE INPUT section. 
It should be noted that the value of the scaling ratio, S “ 2.9, for the 
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5 cm thruotor wao rootrictod by computer CPU time aUocationSi The 
modeling of the upotroam region for a 5 cm thruster would be better if 
more stagoa wore used beeauso for 10 stages, S « 3i2 and for 11 stages 
S « 3.5. 

Sample Input 

A typical oot of input parometors for program PLASIM, to simulate 
a 5 cm ion thruster using mercury (Hg) propellant, are, 

Data Card Content 

1 5 cm thruster, non-uniform density distribution, 
basic configuration. 

2 40 150 -2 20 3 9 1 

3 ,025 .60 .08 .50 .05 .7; ' 

A 7.0E+00 3.SF.+00 A.7E-02 6.0E-19 3.34E-23 

5 .75 1,0 0.0 

6 Propagation of a charge-exchange plasma, 5 cm thruster 

7 Distance along beam axis (meters) Radial distance 
(meters) 

A typical set of input parameters for program PLASIM, to simulate 
a 15 cm thruster using mercury (Hg) pi^opcllant, are 
Data Card Content 

1 15 cm thruster, non-uniform density distribution, 
basic configuration. 

2 40 150 -2 20 3 4 1 

3 .075 .60 .12 .30 .63 .85 

4 S.OE+00 2.5E+00 4.7E-02 6.0E-19 3.34E-23 

5 .75 1.0 0.0 

6 Propagation of a charge-exchange plasma, 15 cm thruster. 

Distance along beam axis (meters) Radial distance 
(meters) . 


7 
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A typical sat of input paramaters for program PLASIM to aimuloto a 
30 cm ion thruster using mercury (Hg) propellant arc, 

Data Card Content 

1 30 cm thruster* non-uniform density distribution, 
basic configuration. 

2 40 150 -2 20 3 2 1 

3 .14 .60 .20 .40 2.0 0.9 

4 0.350E+00 0.175E+00 4.7E-02 6.0E-19 3,34E-23 

5 .75 1.0 0.0 

6 Propagation of a charge-exchange plasma, 30 cm 


thruster, 

7 Distance along beam axis (meters) Radial distance 

(meters) . 

The values given on data cards 3 and 4, in ali of the above cases, 
correspond to the experimental values for those quantities. 
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VERIFICATION 

Two scudico, ono analytic and one experlnental, are uaod for 
the verification of the computer code presented herein * The analytic 
study was conducted an a support activity for the development of the 
subject computer code, and is reported in more detail in APPENDIX A. 

Analytic Solution 

The analytic solution is for the case of an infinitely long 
cylindrical ion beam with a uniform production of charge-exchange ions 
along the beam. The density and potential variations are restricted 
to be only radial so that an analytic solution could be obtained in 
a straightforward manner. 

A computer solution was obtained using the uniform ion density 
option of the program (KEY°-1) and is shown in Fig. 8, The dots occur 
every 60 iterations. The radial density and velocity from the analytic 
solution are compared to the radial density and velocity from the simu- 
lation, as a function of radial distance, in Figs. 9 and 10, respectively. 
The agreement between the analytic solution and the computer simulation 
for the velocity and the density is excellent. 

Experimental Solution 

8 9 

Experimental surveys of the plasma density * are shown in Figs, 

11 and 12 along with comparisons to the computer code for 5 cm and 15 am 
thruotera, respectively. The operating conditions used in the experi- 
mental studies were duplicated to obtain the isodenslty contours 
generated by the computer simulation. Figure 13 shows the isodensity 
contours generated by the computer simulation for a 30 cm thruster. No 
experimental data is presently available for a 30 cm thruster. The ion 
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Dimensionless radial distance, p = r/r^ 

Fig. 9. Comparison of radial densities calculated using the 
computer code and analytic solution. 


ORIGINAL PAGE IS 
OE EOOR QUALITY 



Dimensionless ion velocity 


1 



Fig. 10. Comparison of radial velocities calculated using the 
computer code and analytic solution, 
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trajectory directions obtained using the simulation are in good agreement 
with the experimental measurements of Carruth and Brady. 

Figures 14, 15 and 16 show typical ion trajectories generated by 
the simulation using the data in the SAMPhE INPUT section for a 5 cm, 

15 cm and 30 cm thruster, respectively. . The dots occur every 60 itera- 
tions. 

Figure 17 is a simulation of the 15 cm thruster for use in comparing 
the computer code contained herein to that discussed in relation to 
Pig. 8 of the October 1980 report.^ 

Limitations in Use 

Major factors affecting the accuracy of the simulation obtained are 
the number of ion paths used, the number of iterations performed and the 
time interval used. When the number of ion paths simulated is iiicreased, 
either more iterations must be used or the time interval size decreased 
through use of the variable TIMEMU, This is necessary to keep the 
spacing between the paths comparable to the distance the ion travels in 
one iteration, this is accomplished using the scaling ratio, S, If care 
is not taken in doing this, path crossings will sometimes occur, especially 
among the paths within one beam radius of the thruster. These path 
crossings result from plasma properties changing so rapidly that the 
error in a path location will exceed the local path spacing. The procedure 
used in carrying out the simulation depends on a "laminar" path structure, 
that is, no intersection of paths. The existence of any crossed paths, 
therefore, invalidates local calculations of densities and the other 
related parameters. 

Other limitations arc imposed by core storage allocations, external 
file space, and the CPU time available in a particular machine. 
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Tho present code takes approximately 20 seconds to compile, 150 seconds 
to execute three stages with NUMION»<iO, NUMIT-150 and ICLMRT-1, occupies 
about 140Kg bytes of central memory and writes approximately 125,000 
characters on an external file for three stages, as above, on a Control 
Data Cyber 171 computer, 

It should be noted that the code could be significantly simplified 
and shortened if it were translated to Fortran 77 (Fortran V). 
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APPENDIX A 
ANALYTIC SOLUTION 

A theoretical benehmork is valuable for verification of the 
computer code developed to model the charge-exchange plasma propagation 
in the vicinity of an operating ion thruster. An analytic solution 
is developed herein for that purpose. 

A cylindrical ion beam is assumed with a length vevy much greater 
than r^^, the beam radius. The current density representing positive 
charge- exchange ion production in the beam is assumed uniform along the 
beam. 

In the region exterior to the beam, three basic physical conditions 
are assumed to hold for the ion population and/or the plasma as a whole. 
The first is continuity of ion current represented by 


• 3” = 0 (Al) 

where 3^ is the ion current density. The barometric equation is also 
used to relate plasma density to local potential V 




(A2) 


where is the potential at the reference density n^ and T^ is 
the electron temperature in the region exterior to the beam. Finally, 
energy conservation for singly charged ions is represented by 


I'vl = (v^^ - 2e(V-V^)/m^)^/2 


(A3) 
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where v is the ion velocity and is the ion mass. As a boundary 
condition at the beam edge, ions are assumed to have acquired the Bolun 
velocity 


''o “ '’b ■ 

B 


1/2 


(A4) 


where T is the electron temperature in the beam. 

®B 

The ion current density is related to the streaming velocity by 


-y 

nev 


(A5) 


For the ansumed symmetry, the velocity is radial and is 


V => v(r)r . 


(A6) 


In cylindrical coordinates, Eq. (Al) can be written with the substitu- 
tion of Eq. (A5) as 

1 9v(r) 3V(r) > ^ 8n(r) 3V(r) „ 

which can be solved for V(r) by eliminating n(r) and v(r) using Eqs. 
(A2) and (A3) . The solution can then be written as 

(1 - e(V-V^)/kT^)^/^exp(-e(V-V^)/kT^) = r/r^ . (A8) 

The density and velocity can be calculated using Eq. (A8) along with 
Eq. (A3) or Eq. (A4) . 
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The solution in terms of potential, density and velocity is 
displayed in dimensionless form in Figs. lA, 2A, and 3A. 




Dimensionless Radial Distance , r/r^j 


3 Fig. lA. Potential as a function of radial distance. 
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Fig. 2A. Density as a function of radial distance. 
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Fig. 3A. Velocity as a function of radial distance. 
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APPENDIX B 

BEAM CURRENT DENSITY PROFILES 

Datorraination o£ the total rate of charge exchongo ion production 
in the beam volume must take into account both the beam current density 
profile of the ion thruster and the spatial distribution of neutral 
propellant atoms in the path of the ion beam, 

For theoretical calculations a simplified ion beam current density 
distribution is usually chosen and axial symmetry assumed, A closed- 
form solution exists only for the very simplest distribution which has 

19 

been used for conservative estimates of charge exchange ion production, 

In this case the current density is given by a Dirac delta function, 

« Jjj6(r) , (Bl) 


where J. is the total beam current and 6(r) is the two-dimensional Dirac 

D 

delta function. This expression essentially places all of the beam 
current on the thruster axis where the neutral density follows a simple 
function, thus allowing a closed-form solution, 

A current density profile that is more accurate for the larger, 
multipole thrusters is a uniform beam current density 

where r^^ is the beam radius, Another profile that approximates the beam 
of some divergent field thrusters is a parabolic profile 

^b “ 2J^(1 - r /r^ )Trr^ , 


(B3) 
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where t io the diotonec from the beam axlBi 

Iho expreooiono given above for beam current density are norroaliged 
ouch that 


2it If 
Q 0 

A Gauooian profile io also used by some workers whore 

“(r/r 

- (J^/vr^^)e ^ , (B5) 

subject to the normalizations 
2tf fio 

Jjt, " / / j|^rdrd0 . (B6) 

The uniform, parabolic and Gaussian profiles can also accommodate 
a defined rate of beam spreading as it propagates. The simulation 
described herein does not include effects from beam spreading. 

Neutral gas leaving the accelerator system is taken to be in free 
molecular flow so that the effective neutral density at an arbitrary 
point such as that shown in Big. IB is proportional to the subtended 
solid angle of the ion optics as viewed from the point (r,0,z). 
Effective neutral densities were calculated numerically for an r-z 
matrix with a resolution of 0.1 rj^. For the near field, where about 
half of the charge exchange takes place, the calculated densities are 
given in dimensionless fom in Table IB. 

The calculated distribution of neutral propellant along with the 
beam current density profile allow a calculation of charge exchange ion 
production rate per unit length as a function of distance from the ion 
optics. 





Table IB. Density of Neutral Propellant Eflfltix, n(r,z)/n 
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6 « 


As a first approximation, the charge exG..ango production rate can 
be calculated in closed form if the Dirac delta function is used for 
the ion beam profile. The neutral density along the beam axis is given 
by 




n 


*^o,rof 

2 



(B7) 


where n^ is defined as the density that would provide the correct 
loss rate of neutral propellant through an opening of the same area as 
the beam. The neutral loss rate is then 


N. 


! 1 ^ 
’ o 


(B8) 


where T^ is the propellant temperature, is the propellant mass, and 
k is Boltzmann’s constant. The charge exchange ion production rate per 
unit length is thus 


NQg(z) “ J^Qn(z)/e , 


(B9) 


for small total production rates. Integrating to obtain the total 
production rate gives 



J, On -r. /2e 
b^ o,re£ b 


(BIO) 


As the assumed ion beam profile becomes less peaked, progressing 
from a Dirac delta function through Gaussian and parabolic functions to 
a uniform distribution, the production rate of charge exchange ions 
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will diwiulsh as more o£ the beom passes through the peripheral regions 
of lower neutral propellant density. The two extremes in production 
rates as a function of distance thus use a delta function and a uniform 
function for the beam profile. The results for these two extremes ore 
slw-m in Fig. 2B. Table ZB gives the calculated total charge exchange 
ion production rates for these two extremes, as well as the Intermediate 
parabolic ion beam profile, The ion beam profile is clearly not a 
dominant parameter for charge exchange ion production. 

The simulation developed to model tl)e charge exchange plasma propa- 
gation can be modified for an arbitrary input for the production rate 
as a function of uistance. 


Charge exchonge ion production rate 


4.0 - 




ib= 




0.4 0.8 1.2 1.6 

Distance from grids, z/r^j 


Fig. 2B. Charge-exchange Ion production rates as a function of distance 
from the grids for extremes of possible beam current density 
distributions. 
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Table 2B» Charge-Exchange Ion Production Rates for Different Beam 
Current Density Profiles. 



Production Rate (J, Qn _r, /2e) 

o,ref b ' 

1.00 

0.97 


0.94 
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APPENDIX C 

GLOSSARY OF VARIABLES FOR PLASIM 

This glossary contains definitions of the variables used in the 
driver program PLASIM, the subroutines BLOCK DATA, READER, INIT, CALC, 
CALGD, BOUND, and WRIT and the function subroutine DS, including all the 
variables in the COMMON blocks, The variables used in the plotting 
subroutines, VRSPLT, VRSPL, LNPLT, and PLOTW are not defined here, 
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BK; 

BMCURJ 

GEXSECt 

C: 

COSPARi 

GOSPER! 

DELTAX! 

DELTAZ: 

DELTLX: 

DELTLZ; 

DELTRX: 

DELTRZ! 

DELTRZ ; 

DN(1,N) 

DNI(I): 

DNL: 

DNOB: 

DNR: 

DS! 

DSPLL ! 
DSPLR! 
DUMMYP, 

DSPLIP! 

F: 

FPARS 

FPARX: 


Boltzman’s constant (J/K). 

Beam current (AMPS) . 

Charge exchange cross section (METERS SQUARED) . 

Time interval divided by mass of ion (TIME/UMSION) . 

Cosine of angle between line parallel to path and horizontal. 

Cosine of angle between line perpendicular to path and 
horizontal. 

Difference between two x coordinates on present path. 
Difference between two z coordinates on present path. 
Difference between two x coordinates on left path. 

Difference between two z coordinates on left path. 

Difference between two x coordinates on right path. 

Difference between two x coordinates on right path. 

Difference between two z coordinates on right path. 

Density on path I at iteration N. 

Initial density between paths (I-l) and (1). 

Density to left side of current path. 

Constant used in calculation of the densities, combination 
of quantities including; BMGUR, GEXSEG, RB, NUMION, UTIL, 
q and VELNEU. 

Density to right side of present path. 

Perpendicular displacement from present path to boundary. 

Displacement to left hand path. 

Displacement to right hand path. 

DUMMY! , DUMMYR: Dummy variables used in the calculation of 

the intercepts. Contain intermediate results. 

t 

Displacement of ion along path. 

Total force acting on ion. 

Force acting parallel to the path. 

Component of FPAR acting in x direction. 
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FPARZ; 

FPER; 

FPERXs 

FPERZ; 

FX; 

FZ; 

l! 

ICLERR; 

ICLPLTj 

ICLWRT: 

IDEF: 

IDXNEW; 
IFLAGl ! 


IFLAG2 ; 


Gomponenj: of FPAR acting in z direction. 

Force acting perpendicular to the path. 

Component of FPER acting in x direction. 

Component of FPER acting in z direction. 

Component of F acting in x direction (FPERX + FPARX) , 

Component of F acting in z direction (FPERZ + FPARZ) . 

Allows do loop index, II, to be passed through COMMON, 
path index, 

Used to determine if code generated error messages, for 
non-fatal errors, are written out. 

" 0, Write messages. 

= 1, Do not write messages. 

Used to determine which (if any) plotting routine is used 
» 1, Use subroutine VRSPLT. 

« 2, Use subroutine LNPLT. 

- 3, Use both VRSPLT and LNPLT. 

“ anything else, No plots. 

Frequency with which results of CALC are output, Write 
statement in CALC called after every ICLWRT number of 
iterations. 

Used to determine when the x and z coordinates of an ion 
trajectory exit point are defined (every other time). 

New index for right or left-most path. 

Used to determine if intersections were found to both the 
left and the right without using linear extrapolation, 

= 0, Linear extrapolation not used. 

= 1, Linear extrapolation used on left. 

= 2, Linear extrapolation used on right. 

« 3, Linear extrapolation used in both cases. 

Error flag, tells whether or not the neighboring paths 
are active, 

= 0, No error, path active, 

> 0, Error exists, value references program statement 
where error condition originated. 

= -1, Path to left is not active, 

= -2, Path to right is not active. 

= -3, Neither path (right or left) is active, 

(Negative values do not indicate an error condition,) 
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1FLAG3: 


IFLAG4*. 

IFVAR; 

IFl: 

IF2s 

II; 

IL; 

IN; 

INFO(K); 

INITDO; 

INITIT; 

lOUT; 

IPAG; 

IPATHS ; 

IR; 

IS AT; 

ISTAT(I) ; 

ITITL(K); 

ITTOTN; 


Used to determine when there is a boundary on the right or 
left (or both) of the current path. 

o 1, No boundary intersected by the perpendicular to the 
current path on either side, 
o 2, Boundary intersected on left. 

*= 3, Boundary intersected on right. 

« 4, Boundary intersected on both the left and right, 

Trajectory one renormalization flag, 

“ 0, Continue iterating as usual. 

" 1, Recalculate position, velocity and density of ion 
on path 1. 

Dummy variable used as an argument in an if statement. 

A format for alpha-numeric data I/O, 

A format for alpha-numeric data I/O, 

Do loop index. 

Index (I) of left-most active path, 

Device code for input unit, used in read statements. 

Array storing Information describing run. 

Initial do loop index for DO 80 N which calls CALC. 

Initial iteration index (stage dependent) . 

Device code for output unit, used in write statements. 

Page number for output. 

Device code for output to external files. 

Index (I) of right-most active path. 

Used as path status holder for reading IPATHS and in 
defining right boundary for high resolution upstream run. 

Status of path I, 

=* 0, Path is active. 

= 1 to NMAX, Path is not active, value is iteration number 
where boundary was intersected. 

= 8888, Path is not active, error condition. 

Array storing graph title and axis labels. 

Total number of iterations to be performed. 


ORIGINAL PAGE IS 
OE POOR QUALITY 



71 


IW: 

J! 

KEJ 


KEY! 


UBI 

LASTDO! 

LASTIT; 

LR! 

N! 

NIP (I)! 

NITP! 

NL: 

NLMl! 

NMAX: 

NR! 

NRMl! 

NSTAG! 

NSTAGE! 


Number of words required to generate 80 clmraeters, 

■ 8 For computer B with a word size of 10. 

■ 14 For computers with a word size of 6. 

(Used in COMMON /lO/ and Routines BLOCK DATA, READER, WRIT, 
VRSPLT and VRSPL) . 

Used to determine when WRITE-435 is executed (CALC) . 

Galling parameter for subroutine WRIT, 
a 1, Output heading, initial information and data. 

“ 2, Output interim status of main variables. 

* 3, Finish of a pass, output results, start of new pass. 

» 4, Create file of path coordinates. 

« 5, Create file of position - density triplets. 

Used to determine type of run, 

o 0, Finish plots and terminate (2 blank cards will do) . 

> 0, High resolution upstream pass, right most boundary is 
defined using KEY'TH path of file IPATHS (used to 
define ZBOUND) . 

« ~1, First pass, uniform density distribution. 

< -1, Regular run; first pass, normal non-uniform density 
distribution. 

Label containing computer coda nama, used for output ■ 

Dummy variable denoting last value of do loop index. 

Final (last) iteration index (path status dependent) . 

Determines which side is being considered (DS) , 

» 1, Looking to the left. 

° 2, Looking to the right. 

(Working) number of iterations on present path, I, total 
number of iterations for the present stage. 

Total number of completed iterations on path I. 

Iteration pass number, used for output. 

(Working) number of iterations on left hand path. 

NL minus one (NL-1) , used for indexing left-hand path. 

NUMPRE divided by four (NUMFRE / 4) (see WRIT(4)). 

(Working) number of iterations on right hand path. 

LR minus one (NR-1) , used for indexing right-hand path. 

Allows do loop index NSTAGE to be passed through COMMON* 

Program stage number. 
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NSTGMUj 

NTOTST; 

NUMION; 

NUMITi 

NUMPRE; 

NUMTRI; 

NUMl: 

NUM2; 

PI; 

PIOV2; 

Qs 

RB: 

RBOUND; 

RB2; 

RB95N; 

RT; 

SINOV2; 

SINPAR: 

SINPER: 

SLOPEL; 

SLOPEP; 

SLOPER; 

SPACER: 

TELIN; 


72 

Used to add 10 to N after first stage, (N stage multiplier), 
NSTAG » 1; NSTGMU - 0 

NSTAG > 1; NSTGMU - 1 

Maximum number of stages to be run, 

Number of ion paths. 

Maximum number of iterations to be performed on any one 
path during any one stage. 

Number of ion paths (NUMION) from run which created file 
IPATHS (see WRIT (4)). 

Number of triplets output to be external file. 

Number of ion paths plus one (NUMION + 1) . 

Two times the number of ion paths plus one (2*NUMI0N + 1) . 
Geometrically defined constant, 3.14X59265... 

PI over (divided by) 2. 

Elementary unit of charge (C) . 

Radius of the beam (METERS) . 

Radial boundary in positive x direction (METERS) . 

Beam radius squared (RB ** 2) . 

95 percent of the beam radius divided by 2 times the number 
of ion paths (RB * .95 / (2 * NUMION)). 

Radius of the thruster (METERS) . 

SINPER over (divided by) 2. 

Sine of angle between line parallel to path and horizontal. 

Sine of angle between line perpendicular to path and 
horizontal. 

Slope of path on left between two "working" points. 

Slope of line perpendicular to current path at endpoint. 
Slope of path on right between two "working" points. 

Initial distance between paths in uniform distribution. 
Temperature of electrons in the ion beam (EV) . 
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mOUT: 

THETAP} 

THRLEN; 

TIME; 

TIWEMU; 

TTHNEUs 

UMSIONJ 

um; 

VGURR: 

VELBET; 


VELBOH: 

VELBTLj 

VELBTR; 


VELKEUf 

VELTOTj 

VELTT2: 

VELTX(I).* 

VELTZ(I); 

VELX: 

VELZt 

VL: 

VPREV: 

VR; 

X: 

XINT; 


T(5inpRrature of olectrona outside the ion beam (EV) . 

Angle between line with slope SLOPEP and horizontal. 

Thruster 3ength (METERE). 

Time interval* defines iteration stop size. 

Time multiplier, used to define the time otep in terms of 
multiple of RBOUND / (VELBOH * NUMIT) . 

Temperature of thermal neutrals in chamber (EV). 

Mass of ions (propellant) (KILOGRAMS). 

Utilization factor of propellant (part of propellant turned 
into ions) . 

Plasma potential at present point on path I. 

Present total velocity between path under consideration and 
neighboring path. 

Bohm velocity. 

Present total velocity between path under consideration and 
path on left. 

Present total velocity between path under consideration and 
path on right. 

Thermal velocity of the neutrals in the chamber. 

Present total velocity of ion. 

Present total velocity of ion along path 2. 

Present total velocity component in x direction. 

Present total velocity component in z direction. 

Velocity contribution for this iteration along x direction. 
Velocity contribution for this iteration along z direction. 
Plasma potential on left side of present path. 

Plasma potential at previous point on path I. 

Plasma potential on right side of present path. 

X coordinate of point to be tested. 

X intersection point on boundary line. 




XINTL? 

XINTR; 

XI0N(I,1){ 

XI0N(I,2)! 

XION(I,N)J 

XI0N(I,N+1) 

XPOSL! 

XPOSR: 

XVELMU? 

n • 

M t 

ZBOUND; 

ZCURR; 

ZINT; 

ZINTLt 
ZINTR; 
ZION(I,l) : 
ZION(I,2) ; 

ZION(I,N): 

ZION(I,N+l) 

ZMAX; 

ZPREV; 


Ik 

X intGroectlon, to lefti of llnea SLOPEP and Sl»0PEE» 

X inteeocction, to eight, of linoo SLOPEP and SLOPER. 

X cooedinato of ion trajectory exit point I. 

First X coordinate of ion of ter icoving ion boom along 
trajectory I. 

Present x position of ion. 

5 Newly ealeulated x position of ion. 

X-position (coordinate) holf way along the perpendicular 
displocement (DSPLL) to the neighboring path on the left. 

X-position (coordinate) half way along the perpendicular 
displacement (DSPLR) to the neighboring path on the right. 

X velocity multiplier, used to define initial x velocity in 
terms of some multiple of the Bohm velocity. 

Z coordinate of point to be tasted. 

Z boundary to right of thruster. 

Current z increment, used to get ion exit points. 

Z Intersection point on boundary line. 

Z intersection, to left, of lines SLOPEP and SLOPEL. 

Z intersection, to right, of lines SLOPEP and SLOPER. 

Z coordinate of ion trajectory exit point I. 

First z coordinate of ion after leaving ion beam along 
trajectory I. 

Present z position of ion. 

; Newly calculated z position of ion. 

Maximum z value on path KEY of IPATHS, defines ZBOUND for 
high resolution upstream run. 

Previous z increment, used to get ion exit points, 

Z velocity multiplier, used to define Initial s velocity in 
terms of some multiple of the Bohm velocity. 
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? i 
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I 
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' APPENDIX D 

' COMPUTER CODE LISTING FOR PLASIM 



n 


1 

5 


10 


15 


20 


25 


30 


35 


«0 


^5 


50 


55 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM PLASIM (INPUT, OUTPUT, TAPE5 • INPUT, TAPE6 * OUTPUT# 

1 NPAPAM, paths# TAP67 ■ PATHS, DPPUG • OUTPUT) 

nftIVPR RHUTINE Ai|M|.iA)AA4rAAA*A 

PROGRAM OESCRIPTTONI PROGRAMMER - WIUIAM OEININCfR# P - 30 ^ PI 
PUCE notes on REVISIONS HEfiEl {INCLUDE DATE# INITIALS# LOCATION 
AND CHANGE HADE **♦ PLEASE *♦'*') 

THIS IS THE DPIVFP ROUTINE FOR THE SIMULATION OF A CHARGE « 
EXCHANGE PLASMA SUPROUNOING A PPOAD REAM ION THRUSTER • THIS ROUTINE 
OEFINES the PROGRAM STRUCTURE AND STAGING GY CALLING VARIOUS 
SUBROUTINES. RFADFRi THE FIRST CALL IS TO SUBROUTINE READER WHICH 
READS IN THE RUN SPECIFICATIONS AND INFOPMATION. INITI THE NEXT 
CALL IS TO SUPROUTINF INIT WHICH INITIALIZES VARIOUS PARAMETFPS AND 
constants# DEPTNFS the ion Fxn points ALONG THE REAM EDGE AND 
PERFORMS THE FIRST ITFRATION, CALC« THE NEXT CALL IS TO 
SURROUTINE CALC WHICH COMPUTES THE ION POSITIONS ALONG THE 
TRAJECTORIES. CAIC IS CALLED NUMIT TIMES. THIS COMPLETES THE 
FIRST STAGE. (NUMIT IS THE MAXIMUM NUMBER OF ITERATIONS TO BE 
PEPFORMFO on ANV GTVFN path during any PARTICULAR STAGE). 

INFORMATION FOR THF PJPST (NUMIT - 10) ITERATIONS IS WRITTEN TO AN 
external file# paths. rpRE MgHOPY IS THEN REINITIALIXFD BY 
TRANSFEPINg THF results of the final ten iterations to the CORE SPACE 
TOP THE initial TEN. THE RESULTS OF THESE TEN ITERATIONS ACT AS A 
BASE FOP ANOTHER SET nP (NUMIT - 10) ITERATIONS, CALC IS CALLED 
ANOTHER (NUMIT - 10) TT**IS TO fUL CokE STORAGE WITH ANOTHER SET OF 
PATH COORDINATES WHICH COMPLETES THE SECOND STAGE. THEN CORE IS 
REINITIALIZED ARAIN AND ANOTHER STAGE IS RUN, THIS IS DONE UNTIL 
THE DESIRED TOTAL NIIMPfr oF ITERATIONS IS PERFORMED, OR UNTIL 
ALL PATHS INTERSFCT BOllND At> lES . FINALLY# THE PLOT INDICATOR, 
«ICLPLT», IS CHFCKED TO DETERMINE WHICH PLOTTING POUTINE TO USE TO 
OUTPUT A PLOT OF THE PATH COORDINATES, IF ANY, 

♦ NOTE *'•'* 

ALWAYS PUT TWO PLANK CARDS AT THE END OF THE DATA DECK TO 
SIGNAL A STOP, THE. WORD SIZE DEPENDENT VARIABLES ARFf 
IFl# IF?, INFO, ITITL AND IW, 

SEE PL^CK DATA FOR INSTRUCTIONS ON WORD SIZE 
OFPENDFNT V4PIA01ES AND CODE GENERATED ERROR MESSAGES. 

VARIABLE DICTIONARY ♦♦♦♦♦♦♦♦ 

DN(K,M) I OFN’STTY ON PATH K, AT ITERATION M, 

ICLPLT I USED TO OETERMINE WHICH (IF ANY) PLOTTING ROUTINE IS 
USFD .1 , USE SUBROUTINE VPSPLT, 

•7 , USE SUBROUTINE LNPLT. 

- 3 # USE BOTH VRSPLT AND LNPLT, 

- ANYTHING ELSE # NO PLOTS. 

INITOO I initial OD LOOP INDEX FOR ‘DO BO N» WHICH CALLS »CALC. 
rSTAT(K) » STATUS OF PATH K, 

• 0 # PATH IS ACTIVE, 

• X TO NMAV , PATH IS NOT ACTIVE# VALUE IS ITERATION 

NUMBER WHERE BOUNDARY WAS CROSSED, 

- RBBB , PATH IS NPT ACTIVE# ERROR CONDITION, 

KEY J USE'' TO DETERMINE TYPE OF RUN# 

• 0 f FINISH plots and terminate, 
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c 
c 

60 C 
C 
c 

C NSTA6 
C 

6!5 C NSTAGE 
C NTOTST 
C NUMION 
C NUMIT 
C 

70 C XinN(KiM) 

C ZION(K,M) 

C 

C END OF PROGRAM OFSCRIPTION 

C 

75 C PROGRAM OECLARATinN STATEMENTS. 

C BLANK rOMHON FOR lARGF ARRAYS# 10 - INPUT-OUTPUT# PARAM - PARAMETERS. 

c 

COMMON ZI0N(42,151) ,xrnN(Al#151) #VELT2a5i»#VFLTV (151)# 

1 NIP(4l)#r>N(41,l'51),0NI{A£)#ISTAT(Al) 

RO COMMON / 10 / IN#rDI)T#INF0aA)#KEY#ICLPLT#ICLMRT#ITITL(2P># 

2 IPATHS#TW#IF1(2)#IF?(A)#ICLEPR 

C0MM(5N / PARAM / N,NUMlON#NUMIT#RB#RBnUND#RT#TELOUT#BMCUR#UTIL» 

3 TELIN#THRLFN#UHSTON#VELPDH#ZBOUND#IL#IR#PI#RK,Q,PB95N# 

A nNnp#CEXSFC#TTHNFU,TIME#TIMEMU,XVELMU#ZVELHU#NSTAG# 

5 NSTGMU#NTOTST# PT0V2 

C 

C INPUT parameters# SOf C IF irATIONS# INFORMATION. INITIALIZE VARIABLES. 

C CHECK KFY TO SIGNAL STOP, 

C 

RO i CALL READER 

INITDO ■ 2 
NSTAG - 1 

c plotting routine VRSPLT plots from core memory 
95 c IF (KEY .EO, 0) CALL VRSPLT 

IF (KEY .EO. 0) FHL VRSPL 
if (KEY .EO. 0) STOP 1 
C 

100 C PERFORM INITIAl I7ATnN. 

c 

CALL INIT 
C 

C BEGIN STAGING. THE DO LOOP INDEX *NSTAGE' GIVES THE STAGE NUMPFR, 

10‘i C 

DO 100 NSTAGE ■ 1# NTOTST 
NSTAG " NSTA6F 
IF (NSTAG ,EQ. 1) GO TO 70 
INITDO • 10 

no C 

C REINITIALIZATION PPOCEOURF. THE RESULTS FOR THE FINAL TEN ITERATIONS 
C APE TRANSFERED TO THE COPE AREA FOR THE FIRST TEN ITERATIONS. THFSF 
C TEN ITERATIONS ARF USFD AS A BASE FOR ANOTHER (NUMIT - 10) ITERATIONS. 
C THE TOTAL ITEPATTHN oeR PATH COUNTER IS REDUCED PY ONE TO ENSURE 


> 0 » H'IGH RESOLUTION UPSTREAM PASS# RIGHT MOST 

BOUNDARY IS DEFINED USING "KEY«TH« PATH OF FILE 
IPATHS (USED TO DEFINE 7B0UND) . 

• -I# FIRST PASS# UNIFORM OISTRIBUTTON. 

< -I# REGULAR RUN# FIRST PASS# NON-UNIFORM DISTRIBUTION, 
I ALLOWS nn LOOP INDEX »N$TA6E» TO BE PASSED THROUGH 
COMMON, 

I PROGRAM STAGF NUMBER. 

I MAXIMUM NUMPFP OF STAGES TO BE PUN. 

I NUMbFR of ION paths TO BE SIMULATED. 

I MAXIMUM NUMBFR OF ITERATIONS TO BE PERFORMED ON ANY ONE 
PATH DURING ANY ONE STAGF. 

I X position of ion on /ATH K at ITERATION H, 

« Z position of ion ON PATH K AT ITERATION M. 
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U 5 C THAT S'ACH STAGP DOG^ THP PftOPGR NUHBER OF nERATlOHS* 

C 

DO 50 H « I, 10 

DO AO K - 1 » NUHION 

IF (ISTAT(K) .NE. 0 ) 60 TO AO 

X 20 MM • M + 141 

XION(K^M) « yiON(K,MKH 
nnN(K,M) « ?ion(h#mm) 

IF (M .FO. 10 } 00 TO 30 

ON(K»H} « RN(K#MM) 

125 00 TO 40 

30 ON(y,K} * 0.0 

40 CONTINUF 

50 CONTINUE 

00 60 I » NiJMION 
130 NIPU) » NPtT) « 1 

60 CONTINUE 

C 

C ITERATE «NUMIT« TIMES TO COMPLETE THE FIRST STAGE. ITERATE 
C »*NUMIT ^ 10 »« TIMFS TO COMPLETE THE SUCCESSIVE STAGES, THIS 
135 C SECTION CALCULATES THE ION POSITIONS ALONG THE PATHS. 

C 

70 CONTINUE 

DO 80 NN - INITDO; NUMJT 

N " NN 

140 CALL WRITH) 

CALL CALC 
90 CONTINUE 

C 

C WRITE INFORHATinW FOR FIRST (NUMIT - 10 ) ITERATIONS IN COPE 
145 C TO EXTERNAL FILE. REGiN NFW STAGE IP DESIRED. 

V 

g !ti 1^ 4< « %it> A AAA >}< A A 

C VRSPLT AND LNPLT ^LOT PROM CORE MEMORY 
t r it I VftSPi T 

150 C IF nCLPLT ,GT, 1 ) CALL LNPLT 

f A>^^ 4 ii^tt 4 if 4 . 4 rAAAAAAAAAAAAAAAAAAAAAAAAAAA 
CALL WRIT( 5 ) 

100 CONTINUE 

c 

155 C DETERMINE PLOTTING ROUTINE TO BE USED. IF ANY. AFTER PLOTTING CHECK 
C FOP MOPF INPUT IPPAD MORE DATA IF ANY)* 

C 

IF (ICIPLT .EO. 1 ) CALL VRSPL 
IF (ICLPLT *E 0 , 3 ) CALL VRSPL 
160 60 TO I 

END 


CARD NR. SEVERITY DETAILS DIAGNOSIS OF PROBLEM 

1 ANSI this STATEMFNT TYPE IS NON-ANSI. 
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1 ftLOCK DATA 

C 

C «CnMMnN« DATA ENTRY ROUTINE >!<>»■ >«->»' 

C 

5 C PROOPAH DESCRIPTION? PPOGPAHMER - WILLIAM DEININ6ER, 

C REVISIONS? {INCLUDE OATE» INITIALS AND DESCPIRF CHANGE ** PLEASE 
C 

C THIS PROGRAM LOADS DATA INTO LABLED COMMON STORAGE AT COMPILE TIME 
C THROUGH THE DATA STATEMENTS. 

XO C 

C VARIABLE DICTinNAPY + 

C 

C BK I BOLTZMANis CONSTANT (J/K). 

C IFl I A FORMAT FOR ALPHA-NUMERIC DvATA IZO. 

X5 C IF2 I A FORMAT FOR ALPHA-NUMERIC DATA I/O. 

C IN t DEVICE CODE FOR INPUT UNIT, USED IN READ STATEMENTS. 

C IQUT I DEVICE CODE FOR OUTPUT UNIT, USED IN WRITE STATEMENTS. 

C IPATHS I DEVICE CODF FOR OUTPUT TO EXTERNAL FILES. 

C IW I number of WORDS REQUIRED TO GENERATE PO CHARACTERS, 

20 C "8 FOR COMPUTERS WITH A WORD SIZE OF 10 

C • lA FOP COMPUTERS WITH A WORD SIZE OF h 

C (USED IN «COMMON /I0/» AND ROUTINES »*BLOCK DATA”, ”READER”, 

C rtWRlT”, ”VRSPLT” AND ”VRSPL”) 

C PI I GEOMETRICALLY npFiNED CONSTANT, 3 .1A159265. . . 

?5 C Q « elementary UNIT OF CHARGE (C). 

C 

C END OF PROGRAM DESCRIPTION AND DICTIONARY 

C 

C PROGRAM DECLARATION STATEMENTS. 

30 C BLANK COMMON FOR LARGE ARRAYS, IQ - INPUT-OUTPUT, PARAM - PARAMETERS. 
C 

Common ZI0N(A1,131 ) ,yinN(4l,151>#VILTZ(i51) ,VELTX(i5li, 

I NIP(A1),PN(A1,1MU,0NI(A2),ISTAT(A1) 

COMMON / IQ / IN,rOUT,lNFa(XA),KEY,ICLPLT,ICLWRT,nin(2B), 

35 2 IPATHS,IW,IF1(?),IF2(A) ,ICLERR 

COMMON / PARAM / N,NUMI0N,NUMIT,RB,RB0UND,RT,TELDUT,BMCUR,UTIL, 

3 TE L IN, THR LFN, UM SI ON, VEL BOH, ZBDUN D, U,IR, PI, PK,0,PR9 5N, 

4 0N0B,CEXSFC,TTHNEU,TIMF,TIMEMU,XVELMU,ZVELHU,NSTA6, 

5 NSTGMU,NTnTST, PIOV? 

40 C 

C NOTE? WORD SIZF DFPENPFNT VARIABLES? 

C IFl, IF2, TW, INFO, ITITL 

c IW IS The number of words required to generate bo characters. 

C TO convert to machine using different word SIZE, MODIFY ONLY 

45 c IW, IFl ANII IF? IN first TWO DATA STATEMENTS BELOW. 

C THIRD DATA STATFMrNT FOP I/O BUFFERS, FOURTH DATA STATEMENT FOR 
C VALUES OF CONSTANTS. 

C 

DATA IW, IFKl), TFK?) /fl, 6H(8A10), IH / 

50 DATA IFZtl), IF?.(??, TF2{3), IF2(4) /10H( 8 A10/2( 4A, 4H10 ) ), IH ,1H / 

DATA IN, lOHT, IPATHS /5, 6, 7/ 

DATA 8K, Q, PI / 1 . 3RQ(SF-Z3 , 1.602E-19, 3.141593/ 

C 

C information ON PROGRAM 6FNERATED ERRORS MgsSAGES? 

C WHEN ERROR CALLED, FIRST LINE OF ERROR OUTPUT WILL RE OF THE FORM, 

C «i|<t*>tnk FfjpOR NNN ***^f** 
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WHERE 'NHN» IS ONE QP TMp FQUOWING INTEGERS; 

207 SEE SURROUTINE WRIT 

410 SEE FUNCTION SUBROUTINE 0$ 

412 SEh FUNCTION SUBROUTING OS 

521 SEE SUBROUTINE CALCS' 

5ZZ SHE SU&RfUinNE CAUD 

523 SEE subroutine CAUO 

524 SEE SUBROUTINE CAICO 

525 SEE SUBROUTINE CALCO 

526 SEE SUBROUTINE CALCO 

527 SEE SUBROUTINE CALC 

52B SEE SUBROUTINE CALCO 

529 SEE subroutine CALCO 

530 SEE SUBROUTINE CALC 

610 SEE subroutine BOUND 

6X2 SEE SUBROUTINE BOUND 

AND THE REFERENCED SUPROUTINE IS WHERE THE ERROR IS CALLED FROH, 
'FATAL' ERRORS ARE GENERALLY FATAL TO PARTICULAR PATH ONLY. 

ENO 
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I 

5 


10 


15 


ZQ 


25 


30 


35 


40 


45 


§0 


55 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE PEIDPO 

INPUT {PPAO) ROUTINE *♦♦♦♦«*♦ + ♦♦♦ 

PR06RAH OESCRiPTlONt PROOPAHHER « WlUIAN DEININGER» 7 - 9 - 81« 
REVISIONS! (INCLUne INITIALS AND DESCRIBE CHANGE ** PLEASE 

THIS SUBROUTINE READS IN THE RUN DESCRIPTION AND SPEC IF IC ATIONS > 
BOUNDARY SPECIETCATinNS* PROPELLANT AND PLASHA SPECIFICATIONS# GRAPH 
LABLFS AND READS FJLF IPATHS (IF NECESSARY) TO SUPPLY THE NECESSARY 
information FOR HIGH RESOLUTION UPSTREAM PUNS. 


VARIABLE OlfTIONARY 


BMCUR 

CFXSFC 

iclplt 


iclwpt 


ICLERR 


INFO(K) 

ISAT 

ISTATd! 


niTL(K) 

KEY 


NIP(I) 

NHAK 

NTDTST 

MUMl 

NUMION 

NU>»n 

NU'^PPF 

»R 

RBQUND 

RT 


t BEAM CURRENT (AMPS). 

I CHARDF FVfHAMGE CROSS SECTION (METERS SQUARED). 

I USFD TO OFTERMINE WHICH (IF ANY) PLOTTING ROUTINE IS 

USFD» ■ 1 # USE SUBROUTINE VRSPLT. 

• ? # USE SUBROUTINE LNPLT. 

■ 3 # USE BOTH VRSPLT AND LNPLT. 

• ANYTHING ELSE # NO PLOTS 

I FRFOUFNCY WITH WHICH RESULTS OF CALC ARE OUTPUT. WRITE 

statement in calc called after svery «iclwrt« number 
OF iterations. 

I used to DFTFRMIHE if CODE GENERATED ERROR MESSAGES# FOR 
NON «• FATAL FRRDRS# ARE WRITTEN OUT. 

« Op WRITE MESSAGES. 

« i# DO NOT WRITE MESSAGES. 

« ARRAY STOP IMG INFORMATION DESCRIBING RUN, 

} USFD AS PATH STATUS HOLDER FOR READING IPATHS AND IN 
DEFINING RIGHT BOUNDARY FOR HIGH RESOLUTION UPSTBFAM RUN. 
t STATUS Off PATH I# 

- 0 # PATH IS ACTIVE, 

- 1 TO NMAX # PATH IS NOT ACTIVE# VALUE IS ITERATION 

NUMBER WHERE BOUNDARY WAS CROSSED- 
« BBPB # PATH IS NOT ACTIVE# ERROR CONDITION. 

S ARRAY STORING GRAPH TITLE AND AXIS LABLES. 

S USFD TO OFTERMINE TYPE OF RUN# 

* 0 # FINISH PLOTS AND TERMINATE (2 BLANK CARDS WILL DO), 
> 0 # HIGH RFSOLUTION UPSTREAM PASS# RIGHT MOST BOUNDARY 
IS DEFINED USING KEY'TH PATH OF FILE IPATHS (USED 
TO DEFINE ZBOUND), (SEE NOTE BELOW). 

- -1# FIRST PASS# UNIFORM DISTRIBUTION. 

< -1# REGULAR RUN) FIRST PASS# NORMAL NON-UNIFORM 
DISTRIBUTION, 

t total number of COMPLETED ITERATIONS ON PATH I. 
j NUMPRF DIVIDED BY FOUR (NUMPRE / 4) (SEE WRIT(4)). 
t TOTSL MUMPFR nr STAGES TO BF RUN, 
t NUMPFR OF TON PATHS PLUS ONE (NUMION 1 ) , 
l NUMBFP OF ION PATHS, 

j MAXIMUM NUMBER OF ITFRATIQNS TO BF PREFORMED ON ANY 
ONF PATH DURING ANY ONE STAGF, 

I NUMBER DP TON PATHS (NUMION) FROM RUN WHICH CREATED FILE 
IPATHS rSFF URIT(4)), 

I RADIUS OF THF BEAM (METERS). 

« RADIAL BOUNDARY IN POSITIVE X DIRECTION (METERS). 
t radius Of THE THRUSTER (METERS). 


,Vl 
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C 

C 

60 C 
C 
C 
C 
C 

65 C 
C 

c 

c 

c 

70 C 
C 

c 

c 

c 

75 C 
C 


80 


05 C 
C 

G 

C 

90 C 
C 
C 


95 

C 

C 

C 

100 

C 

C 

C 

105 C 
C 


110 


TEUN I TEHPERATUPE OF ELECTRONS IN THE ION BEAM (EV). 

TELOUT « TEMPFPATUPF OF ELECTRONS OUTSIDE THE ION BEAM (EV). 

THRLEN t THRUSTER LENGTH (METERS). 

TIHEMU I TIME MULTIPLIER^ USED TO DEFINE THE TIME STEP IN TERMS 
OF SOME MULTIBLE OF RBOUND / (VELBOH ♦ NUMIT). 

TTHNEU I TEMPFRATUPF OF THERMAL NEUTRALS IN CHAMBER (EV). 

UMSinSN » HASS OF TONS (PROPELLANT) (KILOGRAMS). 

UTIL « UTILIZATION FACTOR OF PROPELLANT (PART OF PROPELLANT 

TURNED INTO IONS). 

XVELNU » X VELOCITY wuLTIPLlEP» USED TO DEFINE INITIAL VELOCITY 
IN TFPMS of some MULTIBLE OF THE fiOHM VELOCITY. 

ZVELMU I Z VELOCITY MULTIPLIER# USED TO DEFINE INITIAL VELOCITY 
IN TPRMS OF SOME MULTIBLE OF THE BOHM VELOCITY. 

ENO OF PROGPAM DFSCPIPTION AND DICTIONARY 

PROGRAM DECLARATION STATEMENTS. 

BLANK COMMON FOR LARRF ARRAYS# 10 - INPUT-OUTPUT# PARAM ■ PARAMETERS. 

COMMON nON(Al,l51)#XION(4l#l51)#VELTZ(151)#VELTX{151)# 

1 NIP(A1)#DN(A1,151)#DNI( A2)#ISTAT(A1) 

COMMON / ID / IN.I0UT#INF0(1A)#KEY#1CLPLT»ICLWRT#ITITL(Z8)# 

2 IPATHS#IW# IFl(2)#IF2(A)#ICLERR 

common / PAPAM / N#NUMION#NUMIT#PB#RBO(IND#RT#TELOUT#BMCUR»UTIL# 

3 TELIN#THRLEN#IJMSI0N#VELBDH#ZB0UND#IL#IR#PI#BK#0,PB95N# 

A PN0B#CEXSEC#TTHNEU#TIME#TIMEMU#XVFLMU#ZVELMU#NSTAG# 

5 NSTGMU#NTnTST, PT0V2 

READ IN 80 COLUMNS OF INFORMATION DESCRIBING RUN. FIRST CHARACTER 
SHOULD BE A BLANK FOR PRINTER LINE CONTROL. 

READ (IN# IFl) (INFO(K). K - 1# IW) 

READ IN RUN SPECIFICATIONS AND PARAMETERS. 

READ (IN# 12) NUMION, NUMIT# KEY# ICLWRT# ICLPLT# NTOTST# ICLERR 

12 FORMAT (7110) 

IF (KEY ,F0. 0) RFTURN 
NUNl ■ NUMION + 1 

READ IN BOUNDARY SPECIFICATIONS. 

READ (IN# 13) PB» RROUN'J# Rt# THRLEN# BMCUR# UTIL 

13 FORMAT (6F10.5) 

READ IN PROPELLANT AND PLASMA SPECIFICATIONS. TEMPERATURES (FIRST 
THREE VARIABLES LISTED) SHOULD BE INPUT IN ELECTRON-VOLTS. THEY WILL 
BE CONVERTED TO THE DESIRED UNITS FOR CALCULATION BY THE CODE. 

READ (IN# 14) TFLIN# TELOUT# TTHNEU# CEXSEC# UMSION 

14 FORMAT (5F.10.3) 

TELIN - TFLIN ♦ 0 
TELOUT ■ TELOUT ♦ 0 
TTHNEU - TTHNFU ♦ 0 


C 

C READ IN THE TIMF MULTIPLIFR AND THE X AND Z VELOCITY MULTIPLIERS 
C (SEE ABOVE DEFINITIONS). 
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U5 


120 


125 


130 


135 


140 


145 


150 


155 


160 


C 

READ (IN> 15^ TIMPMII, XVEUHU/ ZVELMU 
15 FORMAT (3F10.3) 

C 

C READ TWO CARDS POP GRAPH» flO COLUMNS FOR TITLE AND 40 COLUMNS EACH 
C FOR AXIS LABLES, tWTPn SIZE DEPENDENT AREA) 

C 

IW2 ■ IW t 2 

READ (IN, IF?) (ITTTL(K), K - 1, IW2) 

C 

C TEST TO SEE IF THIS IS A HIGH RESOLUTTCIN UPSTREAM RUN. IF NOT, 

C RETURN TO DRIVER, IF SO, READ FILE IPATHS AND SET UP BnLINOARY. 

C 

IF (KEY) 99, 99# 30 
C 

C REWIND FILE IPATHS SO IT CAN BE READ FOR ANOTHER PASS. THEN READ 
C THE NUMBER OF ION PATHS IN IPATHS TO BE RFAD# THE TOTAL NUMBER OF 
C ION PATHS MAKING UP THE FILE IPATHS (NUMPRE ■ NUMIQN FROM BEFORE)# 

C the status of the PATHS AND THE ARRAY OF PATH COORDINATFS. 

(; note ***** 

Q n ■ ■ ■ n » ■ Ml. ■ ■ ■ K ■ ■■ 

C DUE TO THE OPPINITION OF UMAX (SEE WRIT(4))# KEY MUST RE LESS 

C THAN OP FOUAL TO NMAX FOR A HIGH RESOLUTION UPSTREAM PUN. 

CWO ASSIGNMENT HF PFVTCE CODE SHOULD RF MODIFIFD SO AS NOT TO 

CWD INTFPFER WITH STAGING AND WRIT(5). (SEE WPIT(4)) 

C 

30 REWIND IPATHS 

READ (IPATHS) NMAX# NUHPRF 

READ (IPATHS) (ISTAT(I)# I ■ 1# NMAX) 

DO 50 I • 1# NMAX 
ISAT - ISTAT(I) 

READ (IPATHS) (7IDN(I#NN)# XION(I#NN)# NN • 1# ISAT) 

50 CONTINUE 

C 

C DEFINITION OF THF RIGHT BOUNDARY (USED TO DEFINE ZROUNO). PATH FOR 
C RIGHT BOUNDARY IS DFFJNED BY TRAJECTORY ''KEY” OF FILE IPATHS ALONG 
C WITH THE PATH STATUS AND ITERATIONS ON THE PATH, 

c 

ISAT » ISTAT(KFY) 

DO 70 NN « If ISAT 

7I0N(NIIM1, NN) - 7IDN(KEY# NN) 

XION(NUMlf NN) - XI0N(KFY# NN) 


70 

CONTINUE 



ISTAT(NU'^1 ) 

- ISAT 


NTP(NUMl) 

« ISAT 

99 

CONTINUE 



return 


FND 
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C 

C 

c 

5 C 
C 
C 
C 
C 

10 C 

c 

c 

c 

c 

15 C 
C 

c 

c 

c 

20 C 
C 
C 
C 

c 

25 C 
C 

V 

c 

c 

30 C 
C 
C 
C 

c 

35 C 
C 
C 

c 

c 

40 C 
C 

c 

c 

c 

45 C 

C 

C 

e 

50 c 

c 

c 

c 

c 

55 C 

c 

c 


SUBROUTINE INTT 

INITIAUZATION ROUTINE 

PROGRAH DESCRIPTION! PR0PPAMH6R - WUHAN OEININGER# 7 - 2 - PI 
REVISIONS! (INCLUDE DATE# INTTIALS AND DESCRIBE CHANGE ♦♦ PLEASE *♦) 

THIS SUPROUTINF INITmiZFS THE NECESSARY VARIABLES# CALCULATES 
THE COORDINATES OF THF ION TRAJECTORY EXIT POINTS FROM THE BEAN 
EDGE AND PERFORMS THE FIRSt ITERATION. FIRST THE CONSTANTS ARE 
DEFINED# THEN THF ION EXIT POINTS# THE INITIAL Z AND X POSITIONS 
AND !ZBOUND! ARE CALCULATFO. THEN THE TOTAL VELOCITY COMPONENTS 
AND THE ITERATION NUMBER PER PATH APE INITIALIZED. THE NEXT ION 
POSITIONS ARE CALCULATED AND ALL THE ION PATHS ARE SET TO ACTIVE. 

THE INITIAL DENSITIES ARE CALCULATED AND FINALLY THE RESULTS OF THE 
INITIALIZATION AND FIRST ITERATION ARF OUTPUT. 


VARIABLE DICTIONARY 


ON Id) t 

DNOB I 


IDEF I 

IL t 

IR « 

ISAT I 

ia I M f ^ A f 9 


NIP(U I 
HUMl I 
MUH2 t 
PI0V2 t 
RB2 t 
RB95N I 

SPACER f 

time I 
VEL3ET I 

VELBOH ! 
VELNEU t 
VELTXd) ! 
VELTZd) * 
XIQN( I, l> t 
XI0N(I#2) I 

ZBOUNO t 
ZCURP » 
ZIONdd) t 
ZIONdrZ) t 

ZMAX t 


INITIAL DENSITY BETWEEN PATHS (1*1) AND (I). 

CONSTANT USFO IN CALCULATION OF THE DENSITIES# COMBINATION 
OF OUANTITIFS INCLUDING) BMCUR# CEXSEC# RB# NUHIDN# 

UTIL# 0 AND VELNEU. 

USED TO DETERMINE WHEN THE X AND Z COORDINATES OF AN ION 
trajectory EXIT POINT ARE DEFINED (EVERY OTHER TIME). 

INOFX (I) OF LEFT MOST ACTIVE PATH. 

INDEX (t) OF RIGHT HOST ACTIVE PATH. 

ISTATd) FOR KFY'TH PATH IN UPSTREAM PUN, 

STATUS OF PATH !# 

- 0 # PATH IS ACTIVE. 

- 1 TO NUMTT f PATH IS NOT ACTIVE# VALUE IS ITFPATIDN 

NUMBER WHERE BOUNDARY WAS CROSSED. 

- P8BB # PATH IS NOT ACTIVE# ERROR CONDITION. 

TOTAL NUMBER OF ITERATIONS PERFORMED ON PATH I, 

NUMBER OF ION PATHS PLUS ONE (NUMION + 1). 

TWO times THF NUMBER OF ION PATHS PLUS ONE (2+NUMION + 1). 
PI OVER (DIVIDED BY) 2. 

BEAM RADIUS SQUARED <RB ** 2). 

95 PERCENT OF THE BEAM RADIUS DIVIDED BY 2 TIMES THE 
NUMBER OF ION PATHS (RD ♦ .95 / (2 ♦ NUMION)). 

INITIAL DISTANCE BETWEFN PATHS IN UNIFORM DISTRIBUTION. 
TIHF INTERVAL# DEFINES ITERATION STEP SIZE. 

PRESENT TOTAL VELOCITY BETWEEN PATH UNDER CONSIDERATION 
AND NFIGHBORING PATH. 

BOHM VELOCITY. 

IHERMAL VELOCITY OF THE NEUTRALS IN THE CHAMBER. 

PRESENT TOTAL VELOCITY COMPONENT IN X DIRECTION, 

PRESENT TOTAL VELOCITY COMPONENT IN Z DIRECTION. 

X coordinate of ion TRAJECTORY EXIT POINT I. 

FIRST X COORDINATE OF ION AFTER LEAVING ION BEAM ALONG 
trajectory I. 

Z BOUNDARY TO RIGHT OF THRUSTER, 

PRESENT t INCREMENT, USED TO GET ION EXIT POINTS. 

7 COORDINATE OF ION TRAJECTORY EXIT POINT I. 

FIRST 7 coordinate OF ION AFTFR LEAVING ION BEAM ALONG 
TRAJFCTORY I. 

MAXIMUM 7 VALUF ON PATH **KEY« OF IPATHS# DEFINES ZBDUND 


OF. POOR QUALITY 



S'} 


60 


65 


70 


T5 


60 


65 


90 


95 


100 


105 


110 


C 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


FOR MIPH RF90UITI0N UPSTREAM RUN, 

ZPREV I PREVrnUS i increment# used to get ion exit points* 

END OF PROGRAM DFSCRIPTION AND DICTIONARY 
PROGRAM OECLARATTHN STATEMENTS. 

BLANK common FOR UPGF ARRAYS# 10 » INPUT-OUTPUT# PARAM - PARAMETERS. 


COMMON ZI0N{Al.l5n#XlnN{41#15l)#VeiTm51)#VELTX(151)# 

1 NIP{4U#DN(Alil'51)#DNI(92)#ISTAT(Al) 

COMMON / 10 / IN#irUT#INFn(lA)#KEY#ICLPLT#ICLWRT#ITITL(2e)» 

2 IPATHS»IW#IF1(?)#IF2(4) #ICLERR 

COMMON ! PARAM / N#NUMIQN#NUHIT#RB#RBOUND#RT#TEtOUT# BMCUR»UTR# 

3 TELIN#THRLFN»UMSinN»VELB0H#ZB0UN0#IL#IR#PI»«K#0#R895N# 

A 0N0B»CEXSEr»TTHNBU»TIHE#TIMEMU»XVELHU»7VELMU#NSTAG» 

5 NSTGMU#NTnTST#PlOV? 

DEFINE CONSTANTS INCU'OING BOHM VELOCITY, THERMAL VELOCITY OF THF 

NEUTRALS AND TMB TIMF TNTFRVAL. 


IL 

IP 

NlJMl 

NUM2 

PI0V2 

RB2 

RP95N 

VFLBQH 

VELNEU 

7PREV 

DNOB 

h 

TIME 


1 

NIIMION 
NIIMION + 1 

2 * NIIMION 4 1 
PI / ?,0 

RP ? 

RP * 0,95 / (2,0 FLOAT(NUMION)) 

SORT (TELIN / UMSION) 

SORT (TTHNFU / UMSION) 

0.0 

((BMCIIP ** ?) t CEXSEC 1- (1.0 - UTIL)) / ( FLOAT ( NUMIPN ) 
♦ R6 *■ UTIL * VELNEU ♦ (Q ♦♦ ?) * (PI 2.0)) 

TIMpMU A RBOUND / (VELBOH -I- FLOAKNUMIT ♦ NTOTST)) 


CALCULATION OF THF X ANO 7 COORINATES OF THE ION TRAJECTORY EXIT 
POINTS FROM THE BEAM, GIVES A NDN-LINIFORM DISTRIBUTION OF ION EXIT 
POINTS. (2 * NUMTON 4 11 POINTS ARE CALCULATED# THE EVEN POINTS 
ARE USFD AS ION EXIT POINTS. THE LAST VALUE CALCULATED# 

(2 ♦ NUMION 4 IITH POINT# DEFINES ZPOUND. 


00 100 n • 1# NIJM2 

ZCURR ■ 0.? A (RB95N 4 ZPREV - SORT (ZPREV ♦♦ 2 4 RB2)) - 

7 0,5 * RR? / (RP95N 4 ZPREV - SQRT(7PPEV 4A 2 4 RR?)) 

lOEF - MOD (II# 2) 

IF (IDEF ,E0. 0) GO TO 95 
I - (II 4 1) / 2 

XIONd#!)- RP 
7I0N(I.l)- THRLFN 4 ZdlRR 
95 ZPREV • 7CURP 

100 CONTINUE 

ZROUNO • 7inN(NUMl# 1) 


IF THIS IS A HIGH RESOLUTION UPSTREAM RUN, ZBOUNO MUST PE REDEFINED. 
THE LARGEST ZIONO VALUE ON THE KEY»TH PATH IS USED AS ZBOUND. 

IF (KEY) 120# 120, 105 
105 ZMAX - nON (KFY, 1) 
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115 


120 


125 


130 


135 


140 


145 


150 


155 


160 


165 


170 


ISAT - ISTAT(n 
00 110 H ■ 2# rSAT 

IF (7i0N(KFY» H) ,Le. ZMAX) GO TO 110 
ZMAX • 7I0N(NUH1, H) 

110 CONTINUE 

ZBOUND ■ 7MAX 
C 

C IF KEY EQUALS »1» A UNIFORM DENSITY DISTRIBUTION RUN IS DONE. 

C THE INTERVAL BFTUFEN THRLFN AND ZBOUND IS BROKEN INTO NUHION+l 
C EQUAL INTERVALS. THE ION PATHS START AT THE ENO OF EACH INTERVAL. 
C 

120 IF (KEY .NF. -1) GO TO 125 

SPACER ■ (ZROUNO - THRLEN) / FLOAT(NUMl) 

ZI0N(1»1) - THRLEN 4- SPACER 

DO 122 1 - 2r Nl'MTDN 

ZIONdfl) « 7I0N(I-1; 1) 4 SPACER 
122 CONTINUE 
125 CONTINUE 


C 

C initialization nr the total VELOCITY COMPONENTS AND THE TOTAL NUMBER 
C OF ITERATIONS PFR PATH COUNTER (NlPUH, CALCULATE NEXT ION 
C POSITIONS AND SET ALL ION PATHS ACTIVE. 

C 

DO 130 I - 1. NUMION 

VELTX(I) ■ VELP-OM ♦ XVELMU 

VELTZm ■ velroh a ZVELMU 

NlPdl ■ 2 


XiuNd»2) 
ZI0N(d2) 
ISTATd) 
130 CONTINUE 


YlflNdd) + VELTXd) ♦ TIME 
ZIPNdd) + VELTZd) * TIME 
0 


C 

C CALCULATION OF THE INITIAL DENSITIES. 

C THE FOLLOWING GETS THF INITIAL DENSITIES BETWEEN THE PATHS. 
C 


DO 160 I - 1. NUMl 

IF d .EO. 1) GO TO 152 

IF d *F0. NIIMI) go to 156 

VELBET - (SORT (VELTXd-1) ♦♦ 2 + VELTZd-1) ** 2) + 

8 SORT (VELTXd) 2 + VELTZd) ** 2)) / 2.0 

DNI(I)« DNOB / (((XIONd# 1) + XIONd-1# 1)) / 2.0) 

Q * (ZIONfTi 1) - ZI0N(I-1# in t VELBET) 

GO TO 160 

152 DNI(l) ■ DNOB / (XI0N(1# 1) ♦ (7I0N(l, 1) - THRLEN) 

1 ♦ (SORT (VFLTXd) AA 2 + VELTZd) ♦+ 2))) 

GO TO 160 

156 ONKNUMl) - ONOB t (XIONd-1# 1) A (ZI^OUND -ZI0N(I-1, D) 

2 ♦ (SORT (VELTX(I-l) ♦♦ 2 + VELT7(I-1) ♦♦ 2))) 

160 CONTINUE 


C 

C INITIALIZATION OF THF DENSITY ARRAY 
C 

00 100 I ■ Z» NUHl 

0N(I-1# 1) • (ONId) + ONId-D) / 2,0 

180 CONTINUE 
C 

C RETURN TO DRIVFR AFTER OUTPUTTING HEADING# SCHEMATIC OF THRUSTER# 


ORfOINAL PAOE 

OH pool? QUAUTy 


INITIAt PARAMETERS AND THE TWO SETS OF INFORMATION GENERATED BY 
THIS SUBROUTINE. 


TELIN • TELIN t t) 
TEtOUT ■ TELOUT / 0 
TTHNEU • TTHNEU / 0 
CAU WRIT(I) 

TEIIH • TEIIN ♦ 0 
TELOUT m TELOUT * 0 
TTHNEU • TTHNRU ♦ 0 
N • I, 

CALL WRIT(2» 

N V 3 

CALL WRIT(2) 

RETURN 

END 
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C 

C 

C 

5 C 
C 
C 
C 
C 

10 C 
C 
C 
C 
C 

IS C 
C 
C 
C 
C 

zo c 
c 
c 
c 
c 

25 C 
C 
C 
C 

c 

30 C 
C 
C 

c 

c 

35 C 

c 

c 

c 

c 

AO c 
c 
c 
c 
c 

A3 C 
C 
C 
C 
C 

50 C 
C 

c 

c 

c 

55 C 

e 

c 


SUBRnUTIKE CALC 

^fff^U***** CALCULATION ROUTINE 


PROGRAM DESCRIPTION} PROORAHHER " WILLIAM OEININGER/' 5 " 26 - 81 
REVISIONS! (INCLUDE DATE# INITIALS AND DESCRI8E CHANGE ♦♦ PLEASE **) 

THIS SURRDUTINE USES THE ARRAYS IN BLANK COMMON 
ALONG WITH SUBPOUTINF CALCD (GETS DISPLACEMENTS TO RIGHT AND LEFT! 

TO DETERMINE THF NEXT POSITION OF THE ION BEING CONSIDERED. THE 
DENSITIFS AND POTFNTIALS TO THE BIGHT AND LEFT OF THE CURRENT PATH 
ARE CALCULATED FIRST, THEN THE FORCE ACTING PERPENDICULAR TO THF 
CURRENT PATH IS CALCULATED. THE POTENTIALS AND FORCES ACTING 
PARALLEL TO THE CURRENT PATH APE OBTAINED NEXT, THE TOTAL FORCE 
(SUM OF PERPENDICULAR AND PARALLEL COMPONENTS) ACTING ON THE ION IS 
CALCULATED AND THFN THF VPLOCITY COMPONENTS ALONG THE X AND Z AXES 
ARF OBTAINED, FINALLY THF NEXT ION POSITION IS CALCULATED. SUB- 
ROUTINE BOUND IS CALLED TO MAKE SURE THF NEW ION POSITION IS INSIDE 
THE BOUNDARIES, THE RESULTS ARE PRINTED (EVERY «ICLWRT" TIMES) AND 
FLAGl IS CHECKED (SEE CALCD) TO SEE IF INTERSECTIONS WERE FOUND TO 
BOTH THF, RIGHT AND THE LEFT. IF INTERSECTIONS WERE FOUND TO BOTH 
SIDES# PATH I IS ITERATED AGAIN. 


VARIABLE DICTIONARY 


C ! 

COSPAR I 
caSPgK i 

DN(I#N) ! 
DNL I 
DNOB I 

ONR I 
OSPLIP I 
F I 

FPAR I 
FPftn I 
FPA'^7 I 
FPER 1 
FPERX I 
FPERZ I 
FX I 

FZ I 
I I 

IFLAG4 I 


IFVAR t 
J t 
N I 

NIP(l) I 
NSTGHU t 


TIME INTERVAL DIVIDED BY HASS OF ION (TIME f UMStON) 

COSINF OF angle between line parallel TO PATH AND HORIZONTAL 
COSINE n? ANGIE BETWEEN LINE PERPSHOICULAP TO PATH AND 
HORIZONTAL. 


DENSITY ON PATH I AT ITERATION N. 

DENSITY TO LEFT SIDE OF CURRENT PATH. 

CONSTANT U$Fn IN THE DENSITY CALCULATIONS (C IN THE 
REPORTS ) 

DENSITY TO RIGHT SIDE OF CURRENT PATH. 

DISPLACEMENT DF ION ALONG PATH, 

TOTAL FORCE ACTING ON ION. 

FORCE ACTING PARALLEL TO THE PATH, 

COMPONENT OF FPAR ACTING IN X DIRECTION. 

COMPONENT OF FPAR ACTING IN Z DIRECTION. 

FORCE ACTING PERPEND, T^UlAR TO THE PATH. 

COMPONENT OF FPFR ACWNG IN X DIRECTION, 

COMPONENT OF FPER ACTING IN Z DIRECTION. 

COMPONENT OF F ACTING IN X DIRECTION (FPERX ♦ FPARX). 

COMPONENT OF F ACTING IN Z DIRECTION (FPERZ + FPARZ). 

ALLOWS no LOOP INDEX, II, TO BE PASSED THROUGH COMMON, 
PATH INDEX, 

TRAJECTORY ONE RENORMALIZATION FLAG# 

• 0# CONTINUE iterating AS USUAL. 

- I# RFCALCULATE position, VELOCITY AND DENSITY OF ION 
ON PATH 1. 

dummy vapiablf used as one argument in an if statement, 

USED TO determine WHEN WRITE-433 If;. EXECUTED, 

(WORKING) NUMBER OF ITERATIONS ON PRESENT PATH# I# TOTAL 
NUMBER OF ITERATIONS FOR THE PRESENT STAGE. 

TOTAL NUMBFR OF COMPLETED ITERATIONS ON PATH 1. 

USED TO ADD 10 TO 'N» AFTER FIRST STAGE# 

(N STAGF MULTIPLIER) 


f 
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C NSTAG • I t MSTGHU • 0 

C NSTAG > I t NSTGHU » 1 

60 C SXN0V2 I SINPE» OVER (OTVIDEO BY) 2* 

C SINPAR I SIME HP ANGLE BETWEEN LINE PARALLEL TO PATH ANO HORHONTAL. 

C SlNPEft I SING PF ANGLE BETWEEN LING PERPENOICULAR TO PATH AND 

C HORIZONTAL. 

i VCURR I PLaSHA POTGNTTAl AT CURRENT POINT ON PATH I. 

65 C VELBTL I PRESENT TOTAL VELOCITY BETWEEN PATH UNDER CONSIOERATIUN 

C AND PATH ON THE LEPT.» 

C VELBTR I PRESENT TOTAL VELOCITY BETWEEN PATH UNDER CONSIDERATION 

C ANU PATH m THE RIGHT, 

C VELTOT I CURRENT TOTAL VELOCITY OF ION, 

70 C VELTT2 f CURRENT TOTAL VELOCITY OF ION ALONG PATH 2, 

C VELTXm» CURRENT TOTAL VELOCITY COMPONENT IN X DIRECTION. 

C VELTZm» CURRENT TOTAL VELOCITY COMPONENT IN Z DIRECTION* 

C VELX I VELOCITY CONTRI»UTION FOR THIS ITERATION ALONG X DIRECTION. 

C VEL7 I velocity contribution for THIS ITERATION ALONG Z DIRECTION. 

75 C VL I plasma potential ON LEFT SIDE OF CURRENT PATH. 

C VPREV I PLASMA POTFNTTAL AT PPFVIOUS POINT ON PATH I, 

C VR I plasma potential on right SIDE OF CURRENT PATH. 

C XPOSL I X-POSmON (COORDINATE) HALF WAY ALONG THE PERPENDICULAR 

C displacement (DSPLU to the neighboring PATH ON THE LEFT. 

BO C XPOSR I X-POSITION (COORDINATE) HALF WAY ALONG THE PERPENDICULAR 

C DISPLACEMENT (DS^LR) TO THE NEIGHBORING PATH ON THE RIGHT. 

C XION(I*N)t CURP,kht X POSITION OF ION* 

C XI0N(I»N+1)I NEXT (NEW) X POSITION OF ION. 

C ZION(I#N)t CURRENT ? POSITION OF ION* 

05 C ZIONfliN+Xn NPXT (MFW) 7 POSITION OF ION. 

C 

C ♦♦♦ END OF PROGRAM DESPRIPTION AMD DICTIONARY *** 

C 

c program declaration STATEMENTS, 

90 c BLANK COMMON FOR LARGE ARRAYS. 10 » INPUT-OUTPUT# PARAM « PARAMETERS, 
C 

COMMON nON(Al,15l),XinN(Al#15'i)#VELTZ(151).VELTXa51). 

X NIP(AX)fnN(4X#X5X)»0NI(42)#ISTAT(AX) 

COMMON / 10 / IN,IOllT#INFO(lA}#KEY.ICLPLT.ICLWRT#lTnL(28)# 

95 2 IPATHS#IW»IFX(2).IF2(4).ICLERR 

COMMON / PARAM / N#NUMION.NUHn.RB#RBQUNO#RT.TELOUT.BHCUR.UTIL. 

3 TELIN.THRLEN.UMSTON.VELBOHiZBOUND.U# IB/PI#BK#Q#R895N# 

4 dnob.cexsec.tthneu.time.tihemu.xvelmu.zvelmu.nstag# 

5 HSTGMU#NTnTST»PIDV2 

100 c 

C define constants* begin iteration of each PATH, TEST FOR PATH ONE 
C RENORMALIZATION AND SFT THE CURRENT ^WORKING" NUMBER OF ITfBATIONS 
C ON PATH I. HAKE SURE THE TOTAL NUMBER OF ITERATIONS PERFORMED ON 
C PATH I IS NOT TOO LARGE FOR THE CURRENT STAGE ANO SEE IF CURRENT 
105 C PATH IS ACTIVE. 

C HOTGI ANY quantities OPERATED ON BY «AINT»» AND MULTIPLIED 

C OR DIVIDED BY X TIMES SOME POWER OF TEN* ARE BEING 

C truncated to AVOID COMPUTER ROUND-OFF ERROR. 

C 

110 C - AINT (TTMF / (UMSION ♦ l.OE+15)) 

C • C * l.OFflB 

NSTGMU ■ 0 

IF (NSTAG .GT, 1) NSTGMU - X 
IFLAG4 • 0 
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U5 

120 

125 

130 

135 

140 

145 

150 

155 

160 

165 


DO 450 II ■ l» NUHION 
I •• II 

IF (I .PO. 2) IFLAfi4 « 1 

IF (NlPd) .6T. {(NSTAG ♦ NUHIT) - dNSTAG « 1) ♦ 10)H 

1 IFCAG4 • 0 

300 N - NIPm - (NSTAO - 1) ♦ (NUHIT - 10) 

IF (NIP(I) ,GT, {(NSTAfi t NUMIT) - ((NSTAG - 1) ♦ 10)1) 

2 60 TO 450 

IF (ISTAT(I) »HF* 0) 60 TO 450 
C 

C GET THE DISPUCFMFNTS TO THE RIGHT AND LEFT> CHECK THE ERROR 
C FLAG# B0UNOABY INTFR^FCTION FLAG/ HAKE SURE PARTICLE TRAJECTORY 
C PATHS ARE SMOOTH AMO CALCULATE THE NBEOEO TRIGONOMETRIC FUNCTIONS 
C OF THETAP. 

C 

CALL CALCn (I» IFLAGl/ IFLAG2/ IFLAG3/ OSPLR# OSPLL# THETAP) 
IP IIFLAG? ,GT. 0) GO TO 490 

IF (IFLAG3 .LF, 0 .OR. IFLAG3 .GE. 5) GO TO 497 

IF {IFLAGS .F,Q. 4) GO TO 460 

COSPAR • AINT rCOS (PI0V2 - THETAP) ♦ l.OE+03) 

COSPER • AINT (ros (THETAP) * 1.0E+03) 

SINPAP ■ AINT (SIN (PI0V2 - THETAP) * 1.0E^03) 

SINPEP • AINT (SIN (THETAP) ♦ l.OE+03) 

COSPAR •• COSPAP / 1.0E^03 

COSPER - COSPER / l.OE+03 

Sin PAR • SINPAP / i,ogi.ni 

SINPER - SINPEP / 1.0E+03 

C 

C CALCULATION OF THE (EFT anO RIGHT X-POSITIONS HALF WAY ALONG THE 
C PERPENDICULAR OISPLACEHFNTS TO THE NEIGHBORING PATHS 
C (X-COOROINATFS AT CENTER OF DENSITY CELLS). "IF STATEMENT" 

C OETERMINFS THF DIRECTION OF PATH PROPAGATION. 

C 

SIMQV2 ■ SINPER / 2.0 

IF (ZlONfT.N) - 7I0N(I/N.-1)) 314/ 310/ 318 

C 

310 yPOSL • XION(I/N) 

XPOSP ■ YlONfl.N) 

GO TO 320 
C 

314 XPOSL • YION(I.N) - OSPLL ♦ SINOVa 

XPOSR ■ XIONd/N) + DSPLP ♦ SIN0V2 
GO TO 3P0 
C 

318 XPOSL • XIONd.N) + OSPLL * SIN0V2 

XPDSR - XrONd/N) - OSPLR * SIN0V2 
C 

C CALCULATION OF THE TOTAL VELOCITIES BETWEEN THE PRESENT PATH 
C AND the paths TO THF RIGHT AND LEFT. 

C 

320 IF d .PQ. 1) GO TO 322 

VELRTL - (SORT (VELTX(I-i) ♦♦ 2 + VELTZ(I-l) 2) 

3 + SORT (VELTXm ♦♦ 2 + VELTT(I) ** 2)) / 2,0 
IF (T ,F0, NUMIQN) 60 TO 323 

321 VELBTR • (SORT (VFLTXd) ** Z * VELTZ(I) ♦* 2) + 

4 SORT(VELTXd + l) ♦♦ 2 + VELTZ(Ifl) 2)) / 2.0 
GO I'D 32 5 
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3?.2 VELBTL ■ SORT (VELTXd) ♦♦ 2 + VELTZ(I) 2) 

GO TO 3?1 

323 VELRTP - SORT (VECTX(I) ** Z f VELTZdJ ♦♦ 2) 

XT5 C 

C CALOUUTION OP THE OENSTTIES TO THE RIGHT AND IEFT» AND THE 
C AVERAGE DENSITY AT THE CURRENT POINT. HAKE SURE ONL AND DNR 
C ARE NOT ZERO. 

C 

180 325 ONL « AINT (DNOB / {XPOSL ♦ OSPLL * VELRTU) 

ONR • AINT (f)NOB / (XPOSR * OSPLR ♦ VELATR)) 

DN(I#N) » (ONL + DNR) / 2.0 

IF (ONL .LE. 0.0 .OR. DNR .LE. 0.0) GO TO 495 
C 

185 C WHEN A BOUNDARY IS INTERSECTED ON THE LEFT OR RIGHT» THE DISPLACE - 
C HENTS HAVE TO BF CHECKED TO MAKE SURE NO BOUNDARY REPULSION EXISTS. 
C IF THE D2SPLACFHENT FROM THE CURRENT PATH TO THE BOUNDARY IS LESS 
C THEN THF DISPLACEMENT BETWEEN THE CURRENT PATH AND THE NEIGHBORING 
C PATH^ THE PERPENDICULAR FORCE IS ZEROED# OTHERWISE THE PATH 
190 C PROPAGATES AS USUAL. 

C 

IF (IFLAG3 .EO. 1) GO TO 340 



IF (VELT7(I) 

.EO, 

. 0.0) 

GO 

TO 340 


IF (IFLAG^ 

.EO, 

. 2) 60 

TO 

32B 


IF (IFLAG3 

.EO, 

.3) GO 

TO 

335 


GO TO 497 





328 

IF (DSPLL ♦ 

LE. 

OSPLR) 

GO 

TO 330 


GO TP 340 





330 

FPERX • 0,0 






FPERZ - 0,0 






60 TO 370 





335 

IF (DSPLR . 

LF, 

DSPLL) 

60 

TO 330 


C CALCULATION OF THE POTENTIALS TO THE RIGHT AND LEFT. THE FORCE# 
205 C AND ITS CONPONENTS# ACTING PERPENDICULAR TO THE PATH IS THEN 
C CALCULATED USING THE SMALLER OF THE TWO PERPENDICULAR DIS - 
C PLACEMENTS AS THE DIVISOR. 

C 

340 VL ■ AINT ( (ALOG (DNL / DNI(I)) * TELOUT / 0 ) * 

210 5 l.OE+04) 

VP • AINT {(ALOG (DNR / DNKI+D) * TELOUT / 0) Y 

A l.OE+04) 

VL *» VL / 1 .OE + 04 

VR • VR / 1.0E-fr04 

215 IF (DSPLL .LE. DSPLR) GO TO 345 

FPER ■ ATNT ({0 >♦< (VL « VR) / DSPLR) * l.OE+20) 

GO TO 350 

345 FPER - AINT HO * (VL - VR) / DSPLL) ♦ l.OE + 20) 

350 FPER • FPEP / l.OE+20 

220 FPERX • FPER + SINPER 

FPERZ • FPER ♦ COSPER 

C 

C CALCULATION OF THE POTENTIALS AND THE FORCE AND ITS COMPONENTS 
C ACTING PARALLEL TO THE PATH. 

225 C 

370 VCURR - ALOG (DN(I»N) / DN(I#D) ♦ TELOUT / 0 

VPREV ■ ALOG (0N(I#N-1) / DN(I#D) ♦ TELOUT / 0 

DSPLIP ■ SORT ((XION(I#N) - XION(I#N-D) 2 + (ZI0N(I#N) 
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7 - 7I0N(I,N-1)) ♦♦ Z) 

230 FPAR ■ AINT ((0 (VPREV - VCURR) / DSPLIP) ♦ 1.0E420) 

FPAR • PPAP / l.OE+20 

PPARX « FPAR 4 SINPAR 

FPAR7 ■ FPAR ♦ COSPAR 

C 

235 C CALCULATION OF TOTAL FORCF AND COHPONENTS. 

C 

FX ■ FPARX + FPERX 

FZ • FPAR7 + FPERZ 

F • SORT (FV *4 2 + FZ 44 2) 

2A0 C 

C CALCULATION 0^ VELOCITY COMPONENTS FOR THIS ITERATION# TOTAL 

C VELOCITY COMPONENTS AND TOTAL VELOCITY FOR PATH I. 

C 

VELX - FX 4 C 

2A5 VELZ • F7 4 C 

3<?0 VELTX(I)- VELTXCn 4 VELX 

VELTnn- VFLT7U) 4 VELZ 

VELTOT • SORT (VELTXUJ 44 2 4 VELTZH) 44 ?) 

C 

250 C THIS SECTION MOOTFIFS THE VELOCITY ON THE FIRST PATH SO THAT ITS 
C NORMALIZED WITH RFSPECT TO 1,2 TIMES THE VELOCITY ON THE SECOND 

C PATH IN TERMS OF MAGNITUDE. THE DIRECTION OF THE TRAJECTORY IS 

C LEFT UNCHANGFn, TMIS NORMALIZATION ONLY OCCURS WHEN THE VELOCITY 
C ON PATH 1 IS MORC THAN 20 PRECENT GREATER THEN THE VELOCITY ON 
255 C PATH 2. THE VALUE OF 20 PERCENT IS ARRITARY AND IS BASED ON 

C WHAT 6IVFS THF SMOOTHEST TRAJECTORIES. THIS PREVENTS PATH ONE 

C FROM ACCELERATING TO FAST. 

C 

IF (I .EO. 1 .and. IFLAGA .EQ. 1) GO TO 398 
260 GO TO AlO 

398 VELTT2 • SORT (VELTX(I41) 44 2 4 VELTZ(I41) 44 2) 4 1,2 
IF (VELTOT - VFITT2> AlO# AlO# AOO 
AOO VELTX(I) • (VFLTX(I) 4 VELTT2) / VELTOT 

VELTZ(T) ■ (VFLTZ(I) 4 VELTTP) 7 VELTOT 
265 VELTOT « SORT (VELTXlI) 44 2 + VELTZ(I) 44 2» 

C 

C CALCULATION OF THF NEXT ION POSITION (USES LINEAR APPROXIMATION)# 

C make surf its inside THE BOUNDARIES. 

C 

270 AlO XIQN(I#N41) - VELTX(I) 4 TIME 4 XIONdiN) 

ZI0N(I#N41) « VFLT7(1) 4 TIME 4 ZION(I»N) 

CALL BOUND (7iriN{I,N4l}# XI0N(I#N4l)# I) 

C 

C WRITE THF RESULTS EVERY "ICLWRT" TIMES. IF ISTAT(I) IS NON-ZERO# 
275 C WRITE THE RESULTS. INCREASE NIPd) BY ONE FOR NEXT PASS. 

C IF PATH I reached a BOUNDARY OH THIS ITERATION# SET NIP(I) TO ITS 

C FINAL VALUE. IF TO MANY ITERATIONS HAVE OCCUREO SET ISTAT(I) - N. 

C TEST TC SEE IF PATH I NEEDS TO BE ITERATED AGAIN (TEST IFLAGl) AND 

C CHECK TO SEE TF TTFPATTON LIMIT HAS BEEN REACHED. CHECK IFLAGA 
280 C FOP REN0RMALI7ATI0N OF PATH 1. 

C 

A20 IF (TCLWPT ,LE. 0) GO TO AAO 
J - MOO (NIPd), ICLWRT) 

IF dSTAT(I) ,NF. 0) GO TO A25 
285 IF (J .NT, 0) GO TO AAO 
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295 


300 


305 


310 


315 


320 


325 


330 


335 


340 


425 WRITE (inUT, 435) I, H, NIP(I)# ISTAT(I)^ ZI0NU#NH)^ 

0 XI0Nn,N+l), VEITZU)# VELTXm^ VELTQT» DNdiN) 

^35 FORMAT UV, 13, IX, 3(14, 2X), 6(E13,6, 2X)> 

440 NlPd) ■ NlPd) + 1 

IF (ISTAT(I) .NE. 0) NIP(I) « NIP(l) - 2 
IFVAR • (NSTA6 ♦ NllMIT + 1) - ((NSTA6 « 1) ♦ 9) 

IF (NIP(I) .6T, IFVAR) ISTATd) ■ H 
If (I «E0. 1 ,AND* IFLAG4 ,EQ.l) GO TO 445 
TF (IFLAGl .F.O. 0) GO TO 300 
IF (IFLAG4 .EQ. 0) GO TO 450 
I ■ 1 

IF (ISTATd) ,NE, 0) 60 TO 445 

NIPd) - N2P(I) - 1 
GO TO 300 
445 IFLAG4 • 0 
450 CONTINUE 
NSTGMU • 1 
RETURN 
C 

C IF A BOUNDARY IS INTEPSECTED ON BOTH SIDES OF THE PATH (TO THE RIGHT 
C AND LEFT) ALL FORCES ARB SET EQUAL TO ZERO CAUSING THE X AND 2 
C VELOCITY CONTRIBUTIONS FOR THIS ITERATION TO BE ZERO. THE PATH 
C PROPAGATES LINEARLY, THE DENSITY IS TREATED AS A CONSTANT. 

C 

460 VELX - 0,0 
VELZ ■ 0,0 
DNd,N) - 0N(I,N-1) 

GO TO 390 
C 

C ERROR EXITS AND ERROR CONDITIONS *♦*♦♦♦* 

C 

C - ERROR 527 - TFLAG3 IMPROPERLY DEFINED, FATAL. 

C 

C - ERROR 530 - DNL OR DNR EQUAL ZERO OR LESS THEN ZERO. CAUSES 
C VL OR VR TO BLOW UP, FATAL. 

C 

495 IFLAG2 > 325 

WRITE (IDUT, 530) IPLAGI, IFLAG2, IFLAG3, I# N, DNL, DNR, 

9 OSPLL, DSPLR, XPOSL, XPOSR 

GO TO 499 

497 WRITE (lOUT, 327) IFLAGl, IFLAG2, IFLAG3, I, N 
499 ISTAT( I) ■ 8888 
GO TO 420 

527 FORMAT (/,lIX,23H44t**A ERROR 527 /, /, UX, 

1 33HIFLAG3 IMPROPERLY DEFINED ( FATAL) ,/,llX, 

2 27HCALLED FPOf SUBROUTINE CALC, /,11X,8HIFLAG1 »,I5, 

3 lOH IFLAG2 «iI5,10H IFLAG3 »,I5,5H I -,I5,5H N -,15) 

530 format (/,lIX,23K-«'4*+v4 ERROR 530 llX, 

1 45HDNL OR DNR LESS THEN OR EQUAL TO ZERO ( FATAL ), /,11X, 

2 27HCALLF0 FROM SUBROUTINE CALC,X,11X, 

3 8HIFLAG1 «,I5,10H IFLAG2 -,I5,10H IFLAG3 -,15,58 I ", 

4 15, 5H N -,I5,7H ONL -,E10.3,7H DNR ■ ,EIO .3, T, IIX, 

5 9H OSPLL -,F10.3,9H DSPLR -,E10.3,9H XPOSL -,EI0,3, 

6 9H XPOSR -,E10.3) 

FNO 
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X SUBROUTINE CAI.CD (I, IPtAGl» IFLAGZ^ IFLAG3# OSPLR» OSPU, THFTAP) 

C 

C ***^*** CALCULATION OF THF DISPLACEHENTS TO THE RIGHT AND LEFT ***♦♦♦♦ 
C 

5 C PROGRAM DESCRIPTION t PROGRAHER - WILLIAM DEININGER» 5 - 19 - 81 
C REVISIONS I (INCLUDE OATE» INITIALS AND DESCRIBE CHANGE ♦♦PLFASE**) 

C 

C THIS SUBROUTINE USES THE ARRAYS ZIONO# XIONO AND NIP(I) TO GET 
C THE DISPLACEMENTS FROM THE CURRENT PATH TO THE LEFT AND RIGHT HAND 
10 C PATHS. THE DISPLACEMENT TO THE LEFT IS DENOTED «OSPLl« AND TO THE 
C RIGHT "DSPLR”. FIRST THE PERPENDICULAR IS OBTAINED (SLOPEP)» THEN 
C THE LINE FORMED BY THF FINAL TWO POINTS OF THE PATH ON THE LEFT 
C IS OBTAINED (SLOPEL). AND FINALLY THE X AND 7 INTERSECTIONS# TO 
C THE LEFT# ARE CALCIJLATFD (XINTL# ZINTL). TESTS ARE RUN TO MAKE SURE 
15 C THE INTERSECTION POINTS ARE «GOOD*’ AND TO DETERMINE IF LINEAR 
C EXTRAPOLATION W^S USED. THE DISPLACEMENT TO THE LEFT HAND PATH 
C (DSPLL) IS THEN CALCULATED. LIKEWISE FOR THE RIGHT HAND 
r DISPLACEMENT (DSPLR). FUNCTION (SUBROUTINE) DS IS CALLED TO GET 
C THE DISPLACEMENTS FOR THE SPECIAL CASES# E6l BOUNDARIES ON THE 
20 C RIGHT OP LEFT. 

C 


25 


C 

C 

C 

C 

C 

C 

C 

C 


— -IFLAGl IS USED TO DETERMINE IF INTERSECTIONS WERE FOUND TP BOTH 
THE LEFT AND THE RIGHT WITHOUT USING LINEAR EXTRAPOLATION, 

IFLAGl ■ 0 LINEAR EXTRAPOLATION NOT USED. 

IFLAGl ■ I LINEAR EXTRAPOLATION USED ON LEFT. 

IFLAGl • 2 LINEAR EXTRAPOLATION USED ON RIGHT. 

IFLAGl - 3 LINEAR EXTRAPOLATION USED IN BOTH CASES. 

IF LINEAR EXTRAPOLATION IS NOT USED# THE CURRENT PATH NEEDS TO BE 
ITERATED AGAIN DUE TO ITS SMALLER ACCELERATION. 


30 C 

C IFLAG2 IS THE FRROR FLAG AND TELLS WHETHER OR NOT THE NEIGHBORING 

C PATHS ARE ACTIVE. 

C IFLAG2 ■ 0 

C IFLA62 > 0 

35 C 

C IFLA62 • -1 

C IFLAG2 • -? 

C IFLAG2 - 

C 

AO C — -IFLAG3 IS USED 
C OR left (OR BOTH) 

C IFLAG3 - 1 

C 

C IFLAG3 ■ 2 

A5 C IFLAG3 ■ 3 

C IFLAG3 - A 

C 
C 

C THIS SUBROUTINE RETURNS » IFLAGl# IFLAG2# IFLAG3# DSPLL# DSPIR, 


50 

C 

AND THETAP. 










c 


' VARIABLE 

DICTIONARY 







V 

c 

OELTAXi 

OIFFERFNCF 

BETWEEN 

TWO 

X 

COORDINATES 

ON 

CURRENT 

PATH. 

55 

c 

OELTAZi 

DIFFERENCE 

BETWEEN 

TWO 

z 

COORDINATES 

ON 

CURRENT 

PATH. 


c 

OELTLXi 

DIFFERENCE 

BETWEEN 

TWO 

X 

COORDINATES 

ON 

LEFT 

PATH. 


c 

OELTLZt 

DIFFERENCE 

BETWEEN 

TWO 

z 

COORDINATES 

ON 

LEFT 

PATH. 


NO ERROR. 

FRROR EXISTS# VALUE REFERENCES PROGRAM 
STATEMENT WHERE ERROR CONDITION ORGINATED. 

PATH TO LEFT IS NOT ACTIVE. 

PATH TO RIGHT IS NOT ACTIVE. 

NEITHER PATH (RIGHT OR LEFT) IS ACTIVE. 

TO DFTEPMINE WHEN THERE IS A BOUNDARY ON THE RIGHT 
OF THE CURRENT PATH. 

NO BOUNDARY INTERSECTED BY THE PERPENDICULAR TO 
THF CURRENT PATH ON EITHER SIDE. 

BOUNDARY INTERSECTED ON LEFT. 

BOUNDARY INTERSECTED ON RIGHT. 

BOUNDARY INTERSECTED ON BOTH THE LEFT AND RIGHT 
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95 


] 


). 


}' 


I: 


c 

c 

60 C 
C 
c 
c 
c 

65 C 
C 
C 
C 
C 

70 C 
C 
C 
C 
C 

75 C 
C 
C 
C 
C 

80 C 
C 
C 
C 

c 

85 C 
C 
C 

e 


90 


95 


C 

c 

c 

100 


c 

c 

105 C 
C 

c 


110 


c 


OELTRXI 
DELTRZ* 
DSPLL I 
DSPLR I 
DUHMYP> 

IPLAGl! 
IPLACat 
IFLAGSt 
IFVAR I 
NIP(I)I 
N I 

NL I 
NLMl I 
NR I 
NRHl I 
SLOPEL* 
SLOPEPi 
SLOPERi 
THETAPi 
XINTL I 
XINTR I 
ZINTL t 
ZINTR I 


DIFFERENCE BETWEEN TWO X COORDINATES ON RIGHT PATH. 
DIFFERENCE BETWEEN TWO Z COORDINATES ON RIGHT PATH. 

SEE ABOVE COHHENTS. 

SEE ABOVE COMMENTS. 

DUMMYL» DbMMYP I DUMMY VARIABLES USED IN THE CALCULATION 
OF THE INTERCEPTS. CONTAIN INTERMEDIATE RESULTS. 

SEE ABOVE COMMENTS. 

SEE ABOVE COMMENTS. 

SEE ABOVE COMMENTS. 

DUMMY VARIABLE USED AS AN ARGUMENT IN AN IF JiTATEMENT. 
TOTAL NUMBER OF COMPLETED ITERATIONS ON PATH I, 

(WORKING) NUMBFP OF ITFRATIONS ON PRESENT PATH, TOTAL 
NUMBER OF ITERATIONS FOR PRESENT STAGE. 

(WORKING) number of ITERATIONS ON LEFT HAND PATH. 

NL MINUS 1 (Nl - 1), USED FOP INDEXING LEFT HAND PATH. 
(WORKING) NUMBER OF ITERATIONS ON RIGHT HAND PATH. 

NR MINUS 1 (NR - 1), USED FOR INDEXING RIGHT HAND PATH. 
SLOPE OF path on LEFT BETWEEN TWO ♦•WORKING” POINTS. 

SLOPE OF LINE PERPENDICULAR TO CURRENT PATH AT ENDPOINT. 
SLOPE OF PATH ON RIGHT BETWEEN TWO "WORKING” POINTS. 
angle between LINF with slope "SLOPEP” AND HORIZONTAL. 

X intersection, to left, of lines "SLOPEP” AND "SLOPEL”. 

X INTERSECTION, TO RIGHT, OF LINES "SLOPEP” AND "SLOPEP”. 
Z INTERSECTION, TO LEFT, OF LINES "SLOPEP" AND "SLOPEL". 

Z INTERSFCTION, TO RIGHT» OF LINES "SLOPEP" AND "SLOPER". 


end of PROGRAM DFSCPIPTION AND DICTIONARY 


PROGRAM DECLARATION STATEMENTS. 

BLANK COMMON FOR LARGE ARRAYS, 10 ■ INPUT-OUTPUT> PARAM - PARAMETERS. 

COMMON ZION ( Al, 151 ), X ION ( Al, 151 ),VELTZ( 151 ),VELTX (151) 
1,NIP(A1),DN(A1,151),DNI(A2),ISTAT(AI) 

COHMOAl / 10 / IN,IOUT,INFO(XA),KEY,ICLPLT,ICLWRT,IYITL(28), 

2 IPATHS,IW,IF1(2),IF2(A),ICLERR 

COMMON / PARAM / N »NUMIDN,NUMIT,R0*RBOUNO,RT,TELDUT,BHCUR,UTIL, 

3 TELIN,THRLEN,UMSION, VELB0H,ZB0UND,IL,IR,PI,BK,0,RB95N, 

A DNOB,CEXSfc'C,TTHNEU,TIME,TIHEMU,XVELMU,ZVELMU,NSTAG, 

5 NSTGMU,NTnTST,PI0V2 


INITIALIZE NECESSARY VARIABLES (FLAGS). 


IfLAGl ■ 0 
IFLAG2 ■ 0 
IFLAG3 ■ 0 

CALCULATE SLOPE OF LINE PERPENDICULAR TO THE CURRENT PATH AT END- 
POINT OF CURRENT PATH (NEGATIVE RECIPROCAL OF SLOPE BETWEEN LAST 
TWO POINTS ON CURRENT PATH) AND ANGLE THIS LINE HAKES WITH HORIZONTAL 

DELTAZ - ZI0N(T,N-1) - ZION(I,N) 

DELTAX ■ XION(I,N) - XI0N(I,N-1) 

IF (DELTAX .EO. 0.0) DELTAX ■ l.OE-16 
SLOPEP - DELTA? / DELTAX 
THETAP - ATAN(SLOPEP) 

DUHHYP • X10N(I,N) - SLOPEP ♦ ZION(I,N) 




■■ %1H. 
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115 C 
C 
C 

c 

c 

120 C 
C 

c 

c 

c 

125 C 
C 


130 

C 

C 

C 

C 

135 


C 

140 C 
C 
C 

V 

c 

145 C 


150 


155 

C 

c 

c 

c 

160 C 
C 
C 
C 

165 C 


170 


TWO SUBSECTIONS FOLLOW? THE FIRST SUBSECTION OBTAINS THE INTER- 
SECTION POINT AND DISPLACEMENT TO THE LEFT# THE SECOND SUBSECTION 
OBTAINS THE INTERSECTION POINT AND DISPLACEMENT TO THE RIGHT, 


CALCULATIONS FOR LEFT, 


CHECK TO SEE TF THIS IS THE FIRST ACTIVE PATH AND INITIALIZE 
THE COUNTER FOR THE LEFT HAND PATH. THE FIRST ACTIVE PATH MUST 
BE HANDLED AS A SPECIAL CASE (IL - 1>. CHECK TO SEE If LEFT HAND 
PATH IS ACTIVE. 

IF (I - ID 492# 460# 307 

307 NL • NIP(I-l) - (NSTAG - 1) ♦ (NUMIT - 9) 

NLMl ■ NL - 1 

IF (ISTATU-l)) 470# 311# 4*70 

CALCULATE SLOPE OF PATH ON LEFT BETWEEN FINAL TWO (OR TWO 
"WORKING") POINTS, 

311 DELTLZ ■ 7inN(I-l#NL) - ZION( !-l#NLKl) 

DELTLX ■ XIONtI-l#NL) - XION ( I-1#NLM1) 

IF (DELTLZ .EQ. 0.0) DELTLZ - l.OE-16 
SLOPEL • nELTLX / DELTLZ 

HAKE SURE INTERSECTION! CAN BE FOUND# THE SLOPES OF THE TWO LINES 
CAN NOT BE IHF SAME. STATEMENTS 312 AND THOSE RIGHT BE(,OW# 
BACKSTEP THE LEFT HAND PATH ONE ITERATION# ALLOWING A NEW SLOPE 
TO P.E CALCULATED, IF BACKSTEPPING IS NOT NEEDED# THE 
INTERSECTION POINTS APE CALCULATED. 

IF (SLOPEo - SLOPED 313# 312# 313 

312 IF (IFLAG2 .EO. 312) GO TO 481 
IF (IFLAG2 .EO. 311) GO TO 4G1 
NL - NL - 1 

IF (NL .EQ. 1) GO TO 480 
NLMl - NL - 1 
GO TO 311 

313 OUHMYL ■ XI0N(I-1#NL) - SLOPEL ♦ ZI0N(I-1#NL) 

ZINTL ■ (DUHHYL - DUMMY?) / (SLOPEP - SLOPEL) 

XINTL - SLOPEP * ZINTL + DUMMYP 

TEST TO SEE IF TNTFRSFCTION POINTS ARE GOOD, FIRST FIND DIRECTION 
OF PATH PROPGATION, 320 IMPLIES NEGITIVE Z DIRECTION# 325 IMPLIES 
POSITIVE Z DIRECTION# 330 IMPLIES VERTICAL (UP OR DOWN) DIRECTION. 
THEN TESTS APE RUN TO SEE IF THE INTERSECTIONS ARE "GOOD"# IF 
LINEAR EXTRAPOLATION IS USED (SETS IFLAGD# OR IF BACK STEPPING 
IS NEEDED. THESE TESTS ARE RUN IN ALL CASES. 

IF (ZI0N(T-1»NL) - ZIONd-liNLMD) 320# 330# 325 

320 IF (ZIONf T-1#NL) .LE, ZINTL .AND, ZINTL .LE. 

6 ZinN(I-l#NLMin GO TO 355 

IF (ZINTL ,LT. ZI0N(I-1#ND) GO TO 353 
322 IF (ZINTL .GT. Z ION ( I-l# NLMl ) ) GO TO 312 

IFLAG2 • 322 
60 TO 486 
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C 




325 

IF (ZION( I-l,NLMll .LE. ZINTL .AND. ZINTL .LE 

t 



7 

ZI0N(I-1,NL)J GO TO 355 


175 


327 

IF (ZINTL .GT. ZI0N(I-1,NL)) 60 TO 353 
IF (ZINTL .LT. 7iON(I-l,NLHin GO TO 312 
IFLAG2 ■ 327 





GO TO 4fl6 


IBO 

W 

c 

WHEN STATEMENT 330 IS CALLED WE HAVE TO TEST THE X 

COMPONENTS 


c 

TO 

FIND THE X DIRECTION OF PROPAGATION, THEN SEE IF 

THE 


C 

INTERSECTIONS ARE «GOOD*<, IF LINEAR EXTRAPOLATION IS USED, 


c 

p 

OR 

IF BACK STEPPING IS NEEDED, 


105 

u 

330 

IF (XinNa-l,NU - XI0N(I-1,NLH1J) 340, 404, 

345 


u 

340 

IF (XinN(t-l,NL) .LE. XINTL .AND. XINTL .LE 

t 



0 

XinN(I-X,NLMl)) GO TO 355 

IF (XINTL ,LT. XION(I-l,ND) 60 TO 353 


190 


342 

IF (XINTL .GT. XION(i-l,NLMD) GO TO 312 
IFLAG2 ■ 342 
GO TO 406 




345 

IF (XinNa-l,NLHlJ .LE. XINTL .AND. XINTL , 

LE. 

195 


9 

XI0N(I-1,NL>» GO TO 355 



IF (XINTL .GIT. XI0N(I-1#NU) GO TO 353 


347 IF (VTNTL .IT, XION ( I-1#NLM1 ) ) GO TO 312 

I FLAG?. - 347 
GO TO 4fl6 

200 C 

C STATEHENT 350 RFSETS IFLAG2 IF STATEHENT 404 IS REFERENCED. 

C STATEHENT 353 SETS IFLAGl/ STATEHENT 355 CALCULATES THE 

C DTSPLACEMENT TO THE LEFT. MAKE SURE IFLAG2 IS SET PROPERLY. 
C 


205 

350 

IFLAG2 ■ 0 



GO TO 355 


353 

IFLAGl - IFLAGl 4- 1 


355 

DSPLL - SORT ((XINTL - XION(I,NU ♦* 2 


1 (ZINTL - ZION(I»NH 2) 

210 IF (IFLAG2 .6T. 0) IFLAG2 • IFLAG2 - 312 

C 

C CALCULATIONS FOR RIGHT. 

C — 

C 

215 C CHECK TO SEE IF THtS IS THE LAST ACTIVE PATH AND INITIALIZE THE 

C COUNTER FOR THE RIGHT HAND PATH. THE LAST ACTIVE PATH MUST 

C BE handled as A SPECIAL CASE (IR ■ NUMION). CHECK TO SEE IF THE 

C LEFT HAND PATH IS ACTIVE. 

C 

220 358 IF (IR - T) 49?, 4fl, 359 

359 NR « NII>(I+1) *« (NSTAG - 1) ♦ (NUMIT - 9) 

NRMl • NP - 1 

IF (ISTATd + in 472, 361, 472 
C 

225 C CALCULATE SLOPE OF PATH ON RIGHT BETWEEN FINAL TWO (OR TWO 

C ”WORKING»*) POINTS. 

C 

361 OELTRZ - 7I0N(I+1,NR) - ZION ( I+l, NRMl ) 
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230 


233 


240 


245 


250 


255 


260 


265 


270 


275 


280 


285 


C 

C 

C 

C 

C 

C 

C 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 


c 


OELTRX « - XION 

IP (OEITRZ «K0* 0.0) DELTR2 « l.OE-16 
SLOPBR • OPLTRX / OELTRZ 


HAKE SURE INTPRSECTIONS CAN BE FOUND# THE StOPES OF THE TWO LINES 
CAN HOT BE THE SAHE. STATEMENTS 362 BACKSTEP THE RI6HT HAND PATH 
ONE ITERATION# ALLOWING A NEW SLOPE TO BE CALCULATEOv IF 
BACRSTEPPING IS NOT NFFOED, THE INTERSECTION POINTS ARE 
CALCULATED. 


IF (SLOPEF - SLDPER) 363# 362# 363 

362 IF (IFLAG? .6E. 359 .AND. IFLAG2 .LE. 362) GO TO 4B3 
NR - NR - 1 

IF {NR .EO. X) GO TO 482 
NRHl - NR - 1 
GO TO 36X 

363 OUMHYR • XIONn + l#NR) - SLOPER ♦ ZiaN(I + l#NR) 

ZINTR •* (TUMMYR « DUMHYP ) / (SLOPEP - SLOPER) 

XINTR •• SLOPEP ♦ ZINTR + OUHNYP 

TEST TO SEE IF INTERSECTION POINTS ARE GOOD. FIRST FIND DIRECTION 
OP PATH PROPAGATION# 370 IMPLIES NEGITIVE Z DIRECTION# 375 IMPLIES 
POSITIVE 7 DIRECTION# 380 IMPLIES VERTICAL (UP OR DOWN) DIRECTION. 
AS BEFORE# TFSTS ARE RUN TO SEE IF THE INTERSECTIONS ARE "GOOD»»# 

IF LINEAR EXTRAPOLATION IS USED (SETS IFLAGl)# OR IF BACK STEPPING 
IS NEEDED. THFSF TESTS ARE RUN IN ALL CASES. 

IF (ZI0N(T + 1#NR) - naN(I + l#NRMX)) 370# 380# 375 


^ F V 

2 

372 


, iUn. 7TMTD - I C- 




ZI0N{I+1#NRMX)) GO TO 397 
IF (ZINTR .LT. ZI0N(I+1#NR)) GO TO 395 
IF (ZINTR .GT. 7.I0N(I + l#NRMl)) GO TO 362 
IFLAG2 • 372 
GO TO 4BP 


375 IF (7I0N(H-1»NRM1) .LE. ZINTR .AND. ZINTR *LE. 

3 ZI0N(I+1#NR)) GO TO 397 

IF (ZINTR .GT. nON(I + l#NR)) GO TO 395 
377 IF (ZINTR .LT. ZION { Ul#NRMl ) ) GO TO 362 

IFLAG2 • 377 
GO TO 48R 

WHEN statement 3B0 IS CALLED# WE HAVE TO TEST THE X COMPONENTS 
IN THE SAME MANOR AS WHEN STATEMENT 330 WAS CALLED TO LOOK TO 
THE LEFT. NEED TO TEST FOR «GODOw INTERSECTIONS# LINEAR 
EXTRAPOLATION AND RACK STEPPING. 

380 IF (XIONd + l.NR) - XION( I + X#NRMX) ) 385# 490# 390 


385 IF (XION(I+X#NP) .LE. XINTR .AND. XINTR .LE. 

4 XinN(I+X#NRMX)) GO TO 397 

IF (XINTR .LT. XiaN(I+X»NR)) 60 TO 395 
387 IF (XINTR .GT. XIQN( HX#NRMX) ) GO TO 362 

IFLA62 - 387 
GO TO 4RR 
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^ IF ,L6. XINTR ,AND, XINTR ,LE, 

5 XIDN(T+1»N9)) 00 TO 39? 

IF (XINTR »CT, yiON(I + l»NRn GO TO 395 
i IF (XINTR .IT. XIONU + l»NRHin GO TO 362 

IFLA02 ■ 392 
GO TO APR 

STATFHENT 39A RFSFTS IFLA62 IF STATBHeNT A90 IS REFFRENCEO* 
STATEMENT 395 SETS TPLAGl# STATEHENT 397 CALCULATES THE 
DISPLACEMENT Tn THE RIGHT. MAKE SURE IFLAG2 AND IFLAG3 ARE 
PROPERLY SET. 


IFLAG2 • 0 
GO TO 397 

IFLAGl ■ IFLAOl + 2 

OSPLR - SORT nXINTR - XION(I#NH 


(ZINTR ' 
IF (IFLAG? .GT. 
IF (TFLAG3 .EO. 

GO TO AX8 
IF (IFIAG? ,NF. 
IFLAG3 ■ ^ 

GO TO A50 


. nONd^NU ** 2) 

0) IFLAG2 • IFLAG2 - 362 
?J GO TO AOS 

-2 .AND. IFLAGl #NE. 2) < 


TO A50 


DEFINE IFLAG3. 


THE FOLLOWING SECfTON EXAIMES IFLAGl AND IFLAG2 SO THAT IFLAG3 CAN BE 
DEFINED. IFLAG3 IS USED IN CALC TO DETERMINE WHEN BOUNDARY REPULSION 
EXISTS SO THAT THE NESCESSARY FORCES CAN BE ZEROED. STATEMENT A30i 
BOUNDARIES INTERSECTED ON BOTH SIDES, STATEMENT 4351 BOUNDARY 
INTERSECTED ON RIGHT, STATEMENT 440: BOUNDARY INTERSECTED ON LEFT, 

STATEHENT 445» NO BOUNDARIES INTERSECTED. 


IF (IFLAG2) A20, 445, 496 
IF (IFLAGl) 494, 445, 423 

IF (IFLAG2 .EO. -1 .AND. IFLAGl 
IF (IPLAG2 .FO, -1) GO TO 440 



IF 

(TFLAG? 

,EQ. -2 

.AND. IFLAGl 

325 

IF 

(IFLAG2 

.EQ. -2) 

GO TO 435 


IF 

(IFLAG2 

.NE. -3) 

60 TO 496 


GO TO 445 
60 TO 445 


IF 

(IFLAGl 

.EQ. 

1) 

GO 

TO 

440 

IF 

(IFLAGl 

*EQ. 

2) 

GO 

TO 

435 

IF 

(IFLAGl 

.EQ. 

3) 

GO 

TO 

430 


GO TO 494 

430 TFLAG3 • 4 

GO TO 450 

435 IFLAGT ■ 3 

60 TO 450 

440 IFLAG3 ■ 2 

GO TO 450 

445 IFLA63 ■ 1 

450 IF (DSPIL .GT. 0.25^ GO TO 502 
452 IF (OSPLR .GT. 0.25) GO TO 504 
455 RETURN 

CALCULATIONS FOR (Bf)UNnAPY)/{SP ECIAL CASES, 
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345 

350 

353 

360 

365 

370 

375 

380 

385 

390 

395 


C 

C STATEHENTS 460 ANO 461 CALCULATE THE OISPLACEHENT TD THE LEFT 
C OF the first active path to the 80UN0ARY AND TO THE RISHT OF THE 

C LAST ACTIVE PATH TO THE BOUNDARY# RESPECTIVSLY* WHEN 5TATEHENT 

C 461 IS CALLED IFLAG3 HAS TO BE PROPERLY SET. 

C 

460 OSPLL - OS (nON(I#Nb XION{l#N)i SLOPEP# 0# I) 

IF (OSPLL .FO. I.OE420) 60 TO 496 

IFLA63 • 2 
60 TO 358 
C 

461 OSPLR « D5 (7I0N(IiN)# XiaN(I#K)» SLOPEP# I# I> 

IF (DSPLP .BO. 1.0E420) 60 TO 300 

IF (IFLAGP .NF, *•! .AND. IFLA61 .NE. 1) 60 TO 465 

IFLA63 • 4 
60 TO 466 

465 IFLA03 • 1 

466 RETURN 
C 

C OTHER TESTS AN6 ASSIGNMENTS. 

C 

C DEFINE IFLAG? IN CASES WHERE TESTED PATH IS INACTIVE. 

C 

470 IFLA62 - -1 
GO TO 311 

472 IF (IFLA62) 475# 478# 496 

473 I FLAG? ■ “I 
GO TO 361 

478 IFLAG2 - -2 
GO TO 361 
C 

C BACK-STEPPING LOGIC WHFN MORE THEN TEN BACKSTEPS ARE NF.EOED ON 

C A PARTICULAR PATH. ALLOWS USE OF INFORMATION IN CORE STORAGE 

C NOT YET OVER WRITTEN WITH NEW RESULTS. 

C 

C FOR LEFT-HAND PATH# 

C 

480 IFLAG2 • IFLAG2 4 312 

XFVAR • NIP(I-X) - (NSTA6 - U ♦ (NUMIT - 9) 

IF (IFVAR .GF. 140) GO TO 506 
NLHl ■ 141 
GO TO 311 
C 

481 IFVAR ■ NlP(r-l) - (NSTAO - 1) 4 (NUMIT - 9) 

IF (NLHl .LE. (IFVAR 4 2)) 60 TO 506 

NL • NLHl 
NLHl • NLMl - 1 
GO TO 311 
C 

C FOR RIGHT-HAND PATH# 

C 

482 IFLAG2 • IFLA62 4 362 

IFVAP « Wlpfl+l) - (NSTA6 -1)4 (NUMIT - 9) 

IF (IFVAR .GE. 140) GO TO 508 
NRHl ■ 141 


QUimmi 
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GO TO 361 

C 

AS3 IFVAH » NtP(I*l) - (NSTA6 -* 1) * (NUHIT - 9) 

IF (NRHl ,LE* UFVAR ♦ 2>) GO TO SOB 
NR ■ NRHl 

405 NRHl ■ NRHl - I 

GO TO 361 
C 

C SRROR CONOmONS. 

^ WMH MMM MR MkMI NMWMPP 

41 /:; c 

C ***** STATPHENTS 4f»4 THROUGH 508 ARE VARIOUS ERROR EXITS, ***** 

C THE VALUE OF IFLAG2 DETERMINES IF VALUES OUTPUT ARE FOR RIGHT OR 
C LEFT# SINCE VALUE OF TFLAG2 REFERENCES A PROGRAM STATEHENT. 

C 

415 C 
C 
C 
C 
C 

4Z0 C 
C 
C 
C 
C 

4E5 C 
C 
c 

c 
c 

430 C 
C 
C 
C 
C 

435 C 

c 
c 
c 
c 

440 C 
C 
C 

c 

404 IFLAGZ ■ 404 

445 IF tVELTXfl-l) .FO, 0.0 vANO. VELTZ(I-l) ,EQ. 0.0) 

7GQ TO 350 

WRITE (lOUT# 523) IFLAG2# I, HL# N# nONU-l#NL)# ZION(I-l#NLMl )# 

8 XIDN(I-1#ND# XION{I«l#NLHl)» ZINTL# XINTL# SLOPE?# SLOPEL 
GO TO 510 

450 486 WRITF (IDUT# 522) IFLAG2# I# NL» N# ZION( I-l»NL) # ZION( I-l#NLMl ) , 

9 XiaN(I-l»Nl)» XI0N(I-1»NLHI)# ZINTL# XINTL# SLOPE?# SLOPEL 
GO TO 510 

400 WRITE (lOUT# 522) JFLAG2# I# NR# N# Z ION( I+1»NR) # ZI0N( H-1»NRH1)# 
1 XinNII+l»NR)# XIDN(m»NRHl)# ZINTR# XINTR# SLOPEP# SLOPER 
455 GO TO 510 

490 IFLAG2 “ 490 


* ERROR 521 - STATEMENTS 506# 508 MEAN INDICES NL# NR MERE BACK- 
STEPPED (138 TIMES) UNTIL CORE STORAGE NO LONGER 
contained values THAT WERE CALCULATED DURING THE 
PREVIOUS STAGE. CORRECT COHPARISIONS CAN NOT PE 
HADE PAST THIS POINT. (FATAL). 

« ERROR 522 - STATFHFMTS 486# 488 MEAN THAT THE INTERSECTION POINTS 
CQNSIDEBFD IN STATEMENTS 322# 327# 342# 347# 372# 

377# 307 and 392 DO NOT SATISFY ANY LOGICAL 
CRITERION. UNPHYSICAL INTERSECTIONS (FATAL). 

- ERROR 523 - STATEMENTS 484# 490 INDICATE THAT POINTS N AND N-1 

ARE THE SAME. ION ON LEFT OR RIGHT HAND PATH# 
RESPFCTIVILY# DID NOT HOVE. UNPHYSiCAL UNLESS lON 
HAS ZFRO VELOCITY AND THERE IS NO NET FORCE 
ACTING ON IT (FATAL). 

» ERROR 524 - STATEHENT 492 INDICATES THAT I IS LESS THAN ONE OR 
GREATER THEN NUMIDN. SHOULD NOT OCCUR SINCE I IS 
DO LOOP INDEX (FATAL). 

- ERROR 525 - IFLAGl IMPROPERLY DEFINED (FATAL). 

« ERROR 520^- TFLAG? IMPROPERLY DEFINED (FATAL). 

- ERROR 520 - DSPLL OR OSPLR COULD NOT BE DEFINED (FATAL). 

- ERROR 529 « DSPLL OR DSPLR UNUSUALLY LARGE (NON-FATAL). 
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460 

469 

470 

479 

480 
409 
4Q0 

499 

500 
505 
510 


IF (VELTX(HU «EQ. 0.0 .AN0» VgLmi+l) .EQ. 0,0) 

200 TO 394 

WRITE IIOUT# 523) IFLAC2# li MR# N# ZIOHU+liNR)# ZiaN{I+X#NRMl)# 
3 XIQN(I+1»NR1# XION(IH#NRHl)# ZINTR# XINTRi SLOREP# StOPER 
00 TO 510 

492 WRITE (lOUT# 524) I» ll# IR# NUHION 
GO TO 910 

494 WRITE (lOUT# 929) IFUGl# I# N 
GO TO 910 

496 WRITE UOUT# 526) rFLA62» I# N 
00 TO 910 
498 IFIAG2 « 460 

WRITE nOUT# 928) IFLAG2# I# N 
60 TO 910 


500 IFLA62 • 461 

WRITE nOUT# 920) IPIAG2# I# N 
GO TO 910 

902 IP (ICLFRR «E0, 1) 60 TO 452 

IFLAG2 •* 355 

WRITE tIOUT# 529) IFI.AG2# I# N# NL# OSPLU SLOPEL# SLOPFP# XINTL# 

4 ZINTL# DUHMYL# OU?iH¥P# K10N(I»N)# Zl0NtI#N) 
IFLAG2 " 0 

GO TO 452 

904 IF UCLERR ,EQ« 1) GO TO 910 
IFi.AG2 • 397 

WRITE (IOUT» 929) IFLAG2# I# N» NR# DSPLR# SIDPER# SLOPEP# XINTR# 

5 ?IHTP» nuHKYR. DUMHYP. XIDNIIiN); ZIONa>N) 


IFI.A62 •• 0 

506 WRITE UOUT# 521) IFIAG2# I# NL# N» SLOPEP# Sl.OPEL# ZINTL# 

6 XINTL# DUMHYP# DUHMVL# XION(I-l#NU# XIOHa-l»NLMi)# 

7 ZIONn»l#NL)» nON(i-l#NLMl)# DSPLL# THE7AP 
GO TO 510 

500 WRITE <IDUT# 921) IFLAG2# I# NR# KRHl# N# SLOPEP# SLOPER# ZINTR# 

8 XINTR# DUMHYP# DUNMYR# XION( I+1#NR) # X10N< I+1#HRHD# 

9 zion(i+i»nr)# nnNn+i#NPHi)# dsplp# thetap 
GO to 510 

510 RETURN 


C 

C ERROR CONDITION FORMATS# ERROR NUMBER IS FORMAT NUMRER. 
C 

921 FORMAT t/»UX#23H*t*A«»/t ERROR 521 44*4»4# /# /# IIX# 


1 36HML OR NR RACK STEPPED TO FAR (FATAL)#/? UX# 

2 28HCALLFD FROM SUBROUTINE CALCD#/#11X# ‘ 

3 0HIFLA62 •#t9»9H I -#r9jl3H N(L OR R) -#I5#9H N(L OR # 

4 6HR)H1 •#I9#5H N «#15#10H SLOPEP *#E9.3#12H SLOPE(L OR# 

5 5H R) •#E9.3»/illX»16H ZINT(L DR R) -*#E9.3# /14HXINT(L DR# 

6 5H R) «#E9,3#10H DUMHYP -#E9.3#17H OUMHY(L OR R) «»E9.3» 

7 /#11X.33H YIONd OR +) 1# N(L OR R) ) ■#E9,3# 

6 33HXION(I {- OR ♦) 1# N(L OR P)Ml) -#E9.3# 

9 33H ZlONd (- OR +) 1# N(L OR R)) ■# /» IIX# E9*3# 

1 33HZT0Nd (- OP +) 1# N(L OP R)H1) -#E9.3# 

2 16H OSPL(L OR R) -#E9.3#10H THETAP •#E9.3) 

922 FORMAT ( /#11V#23HAAA4*A ERROR 522 /» /# IIX# 

1 3BHUNPMYSTCAL INTERSECTION POINTS ( FATAL)# /# llX# 

2 2PHCALLED FROM SUBROUTINE CALCD#/#11X# 

3 BHIFLAG2 «#I9,9H I ■#I5#13H N(L OR R) *#I5# 

4 5H N -,I9»33H ZIONd (- OR +) I# N(L OR R)) •# £9,3# /#11X# 
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515 


530 


525 


530 


535 


540 


545 


550 


5 35H2tnNU {- OP +) 1# N{L OR ■#E0.3# 

6 33H VTONCI (- OR ♦) I# N(t OR R)) ■#E9»3i//UX# 

/ 33HXI0N(I (« OR ♦) 1/ NtL OR R)Hl) w»E9.3# 

0 16H 7XNm OR ») »#g9.3#16H XINTtL OR PJ ■#P9.3r/inx# 

9 eHSlOPFP -♦F9*3^7H SI^OPEa OR R) 

523 FORMAT { /aiX# aSH***^** ERROR 523 ******, /thin* 

1 45HinN POSmONS ARE THE SAME# NO MOTION (FATAU#/#UX# 

2 2SHCALLF0 FPQH SUBROUTINE CALCO#/#UX# 

3 6HIFLAG2 -#I5#5H I ■#I5#13H NU OR R) •#I5# 

4 6H N «#I5»33H 2I0NU {- OR I, HU OR F)» -#E9.3#/»UX# 

5 33H2I0N(I C- OR +) I# NR DR R)Hl) *#E9«3# 

6 33H XIONd (- OR +) I# Nil. OR R)) •#E9,3#/#UX# 

7 33HXION{I {« OP +) X# N(t OR RlMl) ••/E9.3# 

0 16H 7INTU DP R» •#E9.3#I6H XXNTU OR R) •#E9,3i /#11X# 

9 8HSL0REP «»F9,3#X7H SIOPEJL OR R> «#E9,3) 

524 FORMAT (/#UX,23H**4**# ERROR 524 ♦+♦*#♦#/# /#X1X, 

1 30HDD LOOP INDEX HOOIFIEO (FATAUW#IXX# 

2 2BHCIUE0 FPOH SUBROUTINE CALCO#/#UX# 

3 SHI -»I5»6H U -#I5#6|I IR -#I5#X0H NUHION "#15) 

525 FORMAT (/#llX#23Ht»**** ERROR 525 *»♦♦**#/#/# UX, 

1 33H7FIA61 IMPROPERLY DEFINED t FATAL) #/# XXX# 

? 2I3HCALLF0 FROM SUBROUTINE CALCD#/#XXX» 

3 OHIFLASX -»I5»5H I »i#I5#5H N *»I5) 

526 FORMAT ( /# XXX # 23H* + t*4* ERROR 526 /#XXX# 

X 33HIFLAC2 IMPROPERLY DEFINED ( FATAL) #/» XXX# 

2 28HCALLFD FROM SUBROUTINE CALCO#/#XXX# 

3 8HlFLAfi2 •»J5#5H I •*#I5»5H N •»I5) 

520 FORMAT (.^# xxx» 23H*'»^»f'# ERROK 52B vv^^t/t /tint 

X 43H0SPLL OP DSPLR COULD NOT BE DEFINED (FATAL) »/# XXX» 

2 20HCAUFD FROM SUBROUTINE CALCD#/#XXX, 

3 0HIFLAG2 I -#I5#5H N "»I5) 

529 FORMAT (/#XXX»23H*’**+4t ERROR 529 *♦*♦**#/#/# XXX# 

X 32H0SPU OR DSPLR LARGE ( NON-FATAL) » /#XXX# 

2 20HCALLEO FRpH SUBROUTINE CALCO#/#XXX# 

3 BHIFLA62 •#I5#5H I *»I5#5H N •»I5#X3H NIL OR R) «# 

4 15#X6H DSPLIL OR R) ■#E9.3#X7H SLOPE(L DP -#E9.3# 

5 /#XXX#0HSLOPEP »>#E9*3#X6H XINT(L OR R) ■»£<>, 3# 

6 X6H 7INT(L OP P) -#E9»3#X7H DUHHYH OR R) -#E9.3# 

7 XOH DilMMYP ■#F«»,3#/»UX#,UHXI0N(I#N) •»E9.3# 

6 X3H 710N(1,N) -,E9.3) 

END 
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C 

20 C 
C 
C 
0 
C 

25 C 
C 
C 
C 
C 

30 C 
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35 
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C 
C 
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C 
c 
c 


50 



SUBROUTINE ROUND t7t V# I) 

BOUNDARY CHECK ROUTINE 

PROGRAM OESCRIPTIONJ PROGRAMMER *» HXUIAM DEJNINGER# X - 3 « B2 
REVISIONS* nNCLDOF OATF» INITIALS AND DESCRIBE CHANGE ♦♦ PLEASE 

THIS SUBROUTINE CHECKS THE POINT UiK) TO SEE IF IT LIES INSIDE 
THE DEFINED BOUNDARIES OF THE SIMULATION, IF (Z#X) DOES LIE INSIDE 
THE DEFINED BOUNDARIES; NO CHANGES ARE MADE AND CONTROL IS RETURNED 
TO SUBROUTINE CALC. IF (7;X) LIES OUTSIDE THE DEFINED BOUNDARIES; 
THE PATH STATUS IS SFT FOUAL TO THE ITERATION NUMBER. IN ADDITION# 
IF tr;X) libs on the first or LASI ACTIVE PATH AND LIES OUTSIDE THE 
BOUNDARIES# THE LFFT-MOST (IL? DR RIGHT-MOST (IR) INDEX IS RESET* 

VARIABLE DICTIONARY 

I I PATH INDFX 

loxNEw » NEW Index fop right or left-host path. 

II I DO LOtrP INDEX. 

IL I INDEX OF LEPT-HOST ACTIVE PATH. 

IR * INDEX OF RIGHT-HOST ACTIVE PATH. 

LASTOQ I DUMMY VARIABLE DENOTING LAST VALUE OF DO LOOP INDEX. 

X * X coordinate of point to be tested. 

I * I COORDINATE OF POINT TO BE TESTED. 


*** END OF PROGRAM OFSCRIPTION AND DICTIONARY 


PROGRAM DECLARATTON STATEMENTS, 

BLANK COMMON FOR LARGE ARRAYS# ID - INPUT-OUTPUT# PARAH ■ PARAMETERS, 


CO 


m vrtu » 


m 9 nki i JL % - i«<! & - %iS4 k-iifm <r k 

fc lun \ *T X 7 ^-ri: r r i aw f i ^ \ A r 9 ^ \ ^ 


X NIP(41}#DN(4l#X5X)#0NitA2)#ISTAT(AX) 

COMMON / 10 / IN#inUT,INF0a4)#KEY#lCLPLT#ICLWRT#imL(28)# 

2 IPATHS#IWfIFl(2)#IF2(4)»ICLERR 

COMMON f PAPAM f N#NUMION#NUMIT#RB#RBOUNO#RT#TELGUT#BHCUR#Um# 

3 TELIN#THRLEN,UMSinN,VELBaH#ZDQUND#IL#IR»PI#BK#Q#RD95N# 

4 0ND8#CEXSEr,TTHNEU#TIME#TIHEMU#XVELMU#ZVELMU#NSTAG, 

5 N$TGHU#NTaTST#PlOV2 


TEST TO SEE IF THE X roOPOINATE IS LESS THEN OR EQUAL TO THE BEAM 
RADIUS# OR GREATER THFN HR EQUAL TO RBOUND. TFST TO SEE IF THE Z 
COORDINATE IS LESS THFN OR EQUAL TO 0 (ZERO! DR GREATER THEN OR 
EQUAL TO 7B0UND. FINALLY, TEST TP SEE IF Z IS LESS THAN OR EQUAL TO 
THE THRUSTER LENGTH AND IF X IS LESS THEN OR EQUAL TO THE THRUSTER 
RADIUS. IF ANY OF THE ABOVE TESTS ARE TRUE# SET ISTATII) • N, 
OTHERWISE RETURN TO SUBROUTINE CALC. 

IF (X .LF. RB» GO TO 50 

IF (X ,GB. RBOUND) GO TO 50 

IF {Z .LS, 0.0) GO TO 50 

IF (Z .GE. 7BDUND) GO TO 50 

IF 17 *LE. THRI.EN .AND. X .LE. RT) GO TO 50 

RETURN 

50 ISTATd) - N 

IF (Z#X) IS ON THF FIRST OR LAST ACTIVE PATH AND LIES OUTSIDE THE 
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115 


C OEHNEO BOUNDARIES# RESET THE COORESPONDING INDEX TO THE INDEX OF 
C THE NEXT ACTIVE PATH. 

C 

C FOR LEFT MOST PATH# 

c 

IP (I - JU 200# eo# 120 

80 CONTINUE 

LASTDO « NUMION « IL 
DO 100 II - 1# LASTDO 
IDXNEW • I + II 

IF (ISTAT(inXNFW)) 210# 105# 100 
100 CONTINUE 

105 IL ■ IDXNEM 

RETURN 
C 

C FOR RIGHT HOST PATH, 

Q IBMniWnwaxRWwmiBMaiaaaiawKaaMi 

c 

120 IF n - IR) 160# 125# 200 

125 CONTINUE 

DO 140 II ■ 1# IR 

IDXNEW ■ I - II 

IF IISTATUDXNFW)) 210# 145# 140 
140 CONTINUE 

145 IR - IDXNEW 

160 RETURN 
C 

C ERROR CONDITIONS. 

C ««««— — — 

c 

C STATEHPWT'5 ?QQ THROUftH CDono cvrTc. 

C -"-.-w 

C - ERROR 610 - I,TL PR IR DEFINED INCORRECTLY# FATAL. 

C - ERROR 612 - ISTATJI) IMPROPERLY DEFINED# FATAL. 

C 

200 WRITE (I0UT#610) I# IL# IR 
iSTATin - enee 
GO TO 250 

210 WRITE {I0UT#612) I# IDXNEW# ISTAKD# I STAT ( IDXNEW ) # N# 

6 NlPtn# NIPUDXNEW) 

ISTAT(I) • P8P8 
250 RETURN ■ 

C 

C ERROR CONDITION FORMATS, ERROR NUMBER IS FORMAT NUMBER. 

C 

610 FORMAT ( /# IIX# 23H»>t'+++* FPROP 610 '^^*****t /p /tlUp 

1 39HI# IL OR IR DcFi^-trD j(NCO<^ftECTLY ( FATAL) #/# IIX# 

2 28HCALLED PROM SUBROUTINE BOUND# /# IIX# 3HI -»I5# 

3 6H IL »#I5#6H IR -#I5) 

612 FORMAT ( /# IIX# aSHA*#**#* ERROR 612 /#11X# 

1 34HISTATn IMPROPERLY DEFINED ( FATAL )»/# IIX# 

2 2RHCALLED FROM SUBROUTINE BDUND#/#HX» 

3 3HI ->I5»10f) IDXNEW -flSflZH ISTAT(I) -#I5# 

4 17H ISTATt IDXNEW) ■#I5#5H N »#I5jlOH NIP(I) ■# 

5 I5#15H NIP(TDXNEW) -#I5) 

END 
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SUBROUTINE WRTTtKU) 

C— WRIT PRINTS INPORHATION ABOUT THE SIMULATION 
C— KE«ll OUTPUT HEADING# INITIAL INFORMATION AND DATA 
C— KE-2» OUTPUT INTERIM STATUS OF MAIN VARIABLES 
C— KE-31 FINISH OF A PASS# RESULTS# START OF NEW PASS 
C— KE«4I CREATE FILE OF PATH COORDINATES 
C--KE-5f CREATE FILE OF POSITION - DENSITY TRIPLETS. 

C 

C BLANK COMMON FOR LARGE ARRAYS 

COMMON nDN(4l#l5l)#XIONUl#l51I#VELTZ{12>l)#VELTX(15l)# 

1 NIP(41I#0N(A1»151)#0NI(42>#ISTAT(41) 

COMMON / to / IN#IOUT#INFO{14)#KEY#ICLPLT#ICLWRT#ITITL{20)# 

2 IPATHS#IW#IF1(2),SF2<4)#ICLERR 

COMMON /PARAH/N,NUMI0N#NUHIT#RB#R80UND#RT#T£L0UT#BMCUR#UTU# 

3 TELIN#THRLEN»UM5ION»VEL8OH»ZftOUND#IL#IR#PI#BK#Q#R095N# 

4 ONaB#CEXSFC#TTHNFU#TIHE#TIHEMU#XVELHU#ZVELMU#NSTAG# 

5 NSTGHU#NT0TST#PIDV2 

DATA IPAG#LAB#NPAS ^0# 6HPLAS IM, 1 / 

C 

C— FORMATS 
C 


10 FORMAT! ////»43H0 THIS RUN MAY BE CHARACTERIZED BY INFO*#// ) 

11 F0PMAT(1H1»//#60X#A4#I3#////// ) 

13 F0RMAT(////XTX,33HP LASMA SIHULATIO N»/^/ 

17X#43HA COMPUTER CODE TO DESCRIBE THE PROPAGATION#// 
17X#43H0F A CHARGE-EXCHANGE PLASMA IN THE VICINITY#// 
17X#40H0F AN ELECTRICALLY PROPELLED SPACECRAFT #/// 
17X,45HWRITTEN BY WILLIAM DEININ6ER AND DALE WINDER#// 
17X»33HF0R THE JET PROPULSION LABORATORY#// 

21X#27H( J p L P. 0. NO. 955322 )#/// 
i7X#4i‘n‘HAP0LD R. KAI 


. I i 


27X»21H0EPARTMENT OF PHYSICS#// 

25X#25HC0L0RAD0 STATE UNIVERSITY#// 

26X#23HF0RT COLLINS# CO 80523 #// 

33X#9HFALL 1981) 

r(lHl#//60X,A6#l3#/////#10X#21HSCHEMATIC OF THRUSTER#/// 
11X#3H- 1#/#UX>3HA I#/#11X#3HJ I#/#11X#3HI I# /# 

11X#3H» i# 5X,6HTHRLEN#/»11X#3HJ «#7X»1HV#/# 
llX#llHt — — — — ,43X#6HZB0UND#/#llX#lHt#6X#4HA I#44X#1HV 
/#9X#ftHRB0UND»3X#4H« : # / IIX# IH I , 6X# IH I # 2X# 45 ( IH- ) # /# 
XlX,lHi,6X#2HRT#5X#lHA#40X#lHt,/#llX#2aHi#6X)2HRB#39X#lHl#/ 
4(11X»1HI#6X#IHI#6X#1H»#40X,1H«#/)#13X,1H+#13(4H- . ),/ 

IIX, 5H(0#0)»50X#lHi, /#5(66X,1H<# /)#21X#45aH-)#/, 
2(21X#lHl»/)»13Xr9(lH-)#/#6(13X#lHl#/ )#///) 

15 FORMAT(18X#?0HINITIAL VALUES OF PARAMETERS#/#/#/# 

5 X# 22 HFpnM SFCONO DATA C AR D# # /# 14X# 6HNUM ION# 5X# 5HNUMIT# 
7X#3HKEY#4X,6HICLWPT#4X,6HICLPLT#4X#6HNTDTST#4X#6HICLERR# /# 
10X»7I10#/#/#/#5X#21HFR0M THIRD DATA C APD# # /# IBX# 2HRB# 4X# 
6HRB0UNn,BX#2HRT#4X#6HTHRLEN#5X#5HBHCUR#6X#4HUTIL#/#10X, 
6F10.3»/#/#/#5X#22HFR0M FOURTH DATA CARD##/#15X#5HTELIN# 
4X#6HTEL0UT,4X.6HTTHNEU,4X#6HCEXSEC#4X#6HUMSI0N»/#10X# 
3F10.3#2P10,3,/W#/#5X#21HFR0M FIFTH DATA CAR0##/#14X# 
6HTIMEMU#4X#frHXVELMU#4X#6HZVELMU, /#l0X#3F10*3,/#/#/#5X» 
22HCALCULATED QUANTITIES# #/#16X# 4HTIME»4X#6HVELB0H» 4X# 
6HZBOUND#/#10X#E10.3#2F10.3) 

“(1H1#12(5H -2- )#A6#I3//10X»27HINTERIM STATUS — ITERATION# 
I4#3H 0F,I4,11H ITERATIONS# 


14 FORMAT! 
2 

3 

4 

5 

6 

7 

8 
9 


1 
2 

3 

4 

5 

6 

7 

8 
9 
1 

21 FORMATI 
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3 /W#3X»IH1»4X^IHN|3X^ 5HISTAT»4X#9H2I0NU+J)#6X#9HXl0N(Ulb 

4 6X#8HVElT7(l)»7X»8HVELTX(n»0X»6X#BX#7HOM(I»N)//) 

22 FQi^HATUXf I3> 1X» 2(U| 2X)^ ^(E13.6f 2X)* 15Xi m*b) 

31 F0RHAT(/W#?7H results OF PASS— ITERATIDN^U^SH 0F»I4#9H JTERATIO# 

2 2HNS^/»3X»IHI#4X»IHM#3X,3HNZP/2X»5HISTAT,4X,9HZI0N(N+1)»6X# 

3 9HXlON«N+l)i6X,BHVELTZa)#7X,8KVELTXab8X;6HVELTOT»0X, 

4 7HDN{I»N)»/J 


35 FORHATT Z»/#5X#39H4*4*4+4444*'» BEGIN WPIT(5) / ) 

36 F0RHAT(/#/#^X#39H4*t**<«4.ti^»*4 END MRIT{5) /) 

C— WHAT KIND OF CALL IS IT 

GOTO (l#2#3/4,5)» KF 
C— 1 1 

C — INITIAL STATE— HEiVDING AND DATA 
1 IPAG»rPA6 4 1 

WRITEnOUT#lULA0»IPA6 


WRITFC TQUT^13) 

IPAG • IPAG + 1 
WRITE(I0UT»14)LA0»1PAG 
IPAG - IPAG f 1 
WRITEtlOUT/in LAB» IPAG 
WRITE(I0UT»101 

WRITE(IOUT,It=l) aNFO{K),K-l,IH) 

WRITEdOUraS) nuhidn, numit# key# iclwrt# iclplt# ntdtst# iclerr# 

1 RB» PBOUND# RT# THRLEN# BHCUR# UTIL# 

2 TFLIN# TFLOUT# TTHNEU# CEXSEC# UMSION# 

3 TINEMIJ# XVELMU# ZVELMU# 

4 time, VELBOH# ZBOUND 


RETURN 

C— 2 2 2 

C— THIS SECTION PRINTS THE INTERIM STATUS AT THE NTH ITERATION 
2 IPAG • IPAG + 1 

WRITEdOUT#?!) LAB,IPAG#N#NUMIT 
DO 28 I-l»NUNinN 

28 WRITEdOUT#22) I# N, ISTATd ) # Z ION (I»N J # X ION (I#N ) # VELTZ (I ) # 
? VFLTXd)#nNd#N> 

RETURN 


C— 3 3 3 

C THIS SECTION PRINTS RESULT OF A PASS AT NTH EXTRAPOLATION 

3 NITP - N + (NSTAG - 1) ♦ (NUHIT + 1) - (10 ♦ (NSTAG - D) 

ITTOTN - (NTHTST * NU«IT + 1) - (NTDTST • 1) ♦ 10 

WRITEdOUT»31 ) NITP# ITTOTN 

RETURN 

C*^" 4 4 4 

C THIS SECTION creates FILE OF PATH COORDINATES 

CWO DEVICE CODES SHOULD RE CHANGED SO AS NOT TO INTERFER WITH WRIT(5J 
4 NMAX-IFlX(f'LnAT(NUMT0N)/4. ) 

REWIND IPATHS 
WRITFdPATHS) NMAX,NUMION 
WRITE d PATHS) (I STATd ) » I- 1 #NMAX ) 

DO 44 I-1#NMAX 
ISAT-ISTATd) 

44 WRITFdPATHS) ( ZION (I »NN ) # X ION (I#NN ) #NN-1# I SAT ) 

RETURN 


C 

C— 5 5 5 

C 

C (WRIT(5) WRITTEN BY WILLIAM DEl'NINGER) 
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C THIS SECTION WRITES INPORHATION FRON THE FIRST INUHIT « 9) 

C ITERATIONS IN CORE HBHORV TO AN EXTERNAL FILE. THE INFORMATION 
C IS STORED ON THE EXTERNAL FILE IN *'TRIPLETS'*» EACH TRIPLET CONTAINS 
C ONE VALUE EACH POR "XIONO"# ♦‘IIONO" AND •’dnO**. THIS IS DONE IN 
C PATH MAJOR ORDFR. IN OTHER WORDS# ALL THE DESIRED RESULTS FOR ONE 
C PATH ARE OUTPUT BEFORE OUTPUTTING ANT RESULTS FOR NEIGHBORING PATHS. 
C FIRST THE PATH STATUS IS CHECKED TO MAKE SURE NO ERROR CONDITIONS 
C WERE SET DURING EXECUTION, THEN THE INITIAL ITERATION INDEX IS SET 

C EQUAL TO ONE FOR THE FIRST STAGE AND 10 FOR ALL STAGES THERE AFTER. 

C THE FINAL ITERATION INDEX IS COMPUTED# FOR WHICH THE MAXIMUM VALUE 
C £S (HUMIT + IJ AND OCCURS IF THE PATH IS STILL ACTIVE. IF THE PATH 

C IS ACTIVE (ISTAT • 0)# THE NUMBER OF TRIPLETS BECOMES (NUMIT - 9», 

C IF ISTAT IS GREATER THEN 7ER0 AND BECAME GREATER THEN ZERO IN 
C THE CURRENT STAGE# THE NUMBER OF TRIPLETS BECOMES (ISTAT - 1). 

C IF ISTAT BECAME NON-ZERO IN A PREVIOUS STAGE WE CONSIDER THE NEXT 
C PATH, 

C AFTER CALCULATING THE NUMBER OF TRIPLETS# THE PATH NUMBER AND 

C NUMBER OF TRIPLETS ARE OUTPUT TO THE EXTERNAL FILE. THEN THE 
C TRIPLETS ARE OUTPUT. THE NEXT PATH IS THEN CONSIDERED# ETC. 

C 

C - ERROR 207 - PATH STATUS IMPROPERLY DEFINED# FATAL, 

C 

8 IF (NSTAG ,E0. 1) RFWINO IPATHS 
WRITE (IOUT#35) 

DO 100 IS - 1# NUMION 

IF (ISTATdS) - BBBfl) 25# 100# 95 
25 INITIT • 1 

IF (NSTAG ,GT, 1) INITIT - 10 

LASTIT - NldlS) - (NSTAG - 1) ♦ (NUMIT - INITIT! 

IF (ISTAT(IS)) 95# 50# 40 

40 XfVaR • (((NSTAG - 1) 4 NUMIT) ♦ 1 - ((NSTAG - 2) ♦ 

1 INITIT) ♦ NSTGMU) 

IF (NIP(IS) .LE. IFVAR) GO TO 100 
LASTIT • ISTATdS) - 1 
GO TO 60 

50 LASTIT - LASTIT - 10 

60 NUMTRI ■ LASTIT 

MDDNTR ■ (NUMTRI / 3) + 1 
write (IPATHS) IS# MODNTR 
DO 90 NS - 1# NUMTRI 



IF (NS ,E0. 1) 

J • MDD (NS# 3) 

GO TO 05 


IP (J .NE. 0) 

GO TO 90 

85 

WRITE (IPATHS) 

XION(ISfNS)# ZI0N(IS#NS)# 

90 

CONTINUE 
GO TO 100 


95 

WRITE (lOUT# 207) 
ISTA' IS) ■ 8889 
RETURN 

IS# ISTAT(IS) 


100 CONTINUE 

WRITE (I0UT»36) 
RETURN 


207 FORMAT (/#11X# 23H«^**^** ERROR 207 ♦♦**♦♦#/# 7# IIX# 

1 38HPATH STATUS IMPROPERLY DEFINED ( FATAL )# /#11X# 

2 29HCALLE0 FROM SUBROUTINE RE ADER# /# IIX# 

3 6HISTAT(#I4#5H ) • #15) 

END 
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50 C 


FUNCTION OS (Z# SUOPlP# LR# II 
♦ BOUNDARY DISPLAiSEMENT ROUTINE 

PROGRAM DESCRIPTIDNf PROGRAMMER - HIUIAM DEININGER# 0-26-01 
REVISIONS! (INCLUDE DATE, INITIALS AND DESCRIBE CHANGE ♦♦ PLEASE 

THIS FUNCTION SUBROUTINE FINDS THE PERPENDICULAR DISPLACEMENT FROM 
THE FIRST OR LAST ACTIVE PATH TO THE BC^yNOARY. CALCD CONSTRUCTS A 
PERPENDICULAR TO THE PATH FROM THE CURRENT POINT (Z#X) MITH SLOPE 
"SLOPEP"* FUNCTION OS THEN EXTRAPOLATES THIS PERPENDICULAR OF SLOPE 
"SLOPEP« TO THE LEFT OR RIGHT (P'^PENDING ON WHETHER VE ARE CONSID - 
ERING THE FIRST OP LAST ACTIVE PATH) AND FINDS THE 7 . AND X INTERCEPTS 
(ZINT AND XINT) ALONG THE BOUNDARY LINE. ZINT AND XINT ARE CHECKED 
TO SEE IF THEY LIE ON OR BETWEEN THF BOUNDARY ENDPOINTS ON THE 
BOUNDARY LINE. TF THEY DO^ THE DISPLACEMENT IS CALCULATED^ IF THEY 
DO NOT» ZINT AND VINT ARF CALCULATED ALONG THE NEXT BOUNDARY LINE 
AND TESTED AGAIN* THIS CONTINUES UNTIL “GOOD** INTERSECTION POINTS 
ARE FOUND OR ALL BOUNDARIES HAVE BEEN CONSIDERED IN WHICH CASE AN 
ERROR MESSAGE IS OUTPUT. THE PERPENDICULAR DISPLACEMENT IS RETURNED 
TO CALCD. 

VARIABLE DTCTIONARY 

DS t PERPENDICULAR DISPLACEMENT FROM CURRENT POINT TO BOUNDARY. 

I ! PATH INDEX. 

LR t DETERMINES WHICH SIDE IS BEING CONSIDERED/ 

■ 1 LOOKING TO THE LEFT. 

- 2 LOOKING TO THE RIGHT. 

X I CURRENT X POSITION ON FIRST OR LAST ACTIVE PATH (XIONd/NI). 

XINT 1 X INTERSECTION POINT ON BOUNDARY LINE. 

Z J CURRENT Z POSTTION ON FIRST OR LAST ACTIVE PATH (ZIDN{I/N)). 

ZINT t Z INTERSECTION POINT IN BOUNDARY LINE. 

END OF PROGRAM DESCRIPTION AND DICTIONARY 

PROGRAM DECLARATION STATEMENTS. 

BLANK COMMON FOR LARGE ARRAYS/ ID ■ INPUT-OUTPUT/ PARAM - PARAMETERS. 

COMMON ZIQN(41/15l)/XI0N{41/151)/VELT7(151)/VELTX{151)/ 

1 NIP(41)/DN(41/151)/DNK42)/ISTAT(41) 

COMMON / 10 / IN,inUT,INF0(14)/KEY/ICLPLT/ICLWRT/lTnL(28)/ 

2 IPATHS/TW/IF1(2)/IF2(4)/ICLERR 

COMMON / PARAM / N/NUMIQN/NUMIT/RB/RBOUND/RT/TELQUT/0MCUR/UTIL/ 

3 TELIN/THRLFN/UMSI0N/VELBOH/ZBOUND/IL/IR/PI/BK/O/RB95N/ 

4 DNOB/CEXSEC/TTHNEU/TIME/TIMEMU/XVEIMU/ZVELMU/NSTAG/ 

5 NSTGMU/NTOTST/ PI0V2 

DETERMINE WHETHER THE DISTANCE ON THE RIGHT OR THE LEFT IS DESIRED. 

IF (LR .EO. 2) GO TO 200 


C 

C CALCULATIONS FOR LEFT. 

55 C 

C TEST iOR INTERSECTIONS ALONG THE END OF THE THRUSTER BETWEEN THE 
C REAM EDGE (THRLEN. RB) AND THE THRUSTER CORNER (THPLEN/ RT). (FIRST 
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C THE INTERSECTinMS ARE FOUND ALONG THE BOUNDARY LINE AND THEN TESTED 
C TO SEE IF THEY LIE ON OR BETWEEN THE BOUNDARY ENDPOINTS ON THE 
60 C BOUNDARY LINE.) (EQUATION OF BOUNDARY LINEI Z » THRLEN) 

C 

100 IF (SLOPFP .LT. i-I.OE+10) GO TO UO 
ZINT • THRLEN 

XIHT • (THRLEN - Z) ♦ SLOPEP + X 
65 IF (XINT .6E. PB .AND. XINT .LE, RT) GO TO 300 

C 

C TEST FOR INTERSFCTIONS ALONG THE EDGE OF THE THRUSTER BETWEEN THE 
C THRUSTER CORNER (THRLEN# RT) AND THE SPACE CRAFT WALL (0# RT). 

C (EQUATION OF BOUNDARY LINEI X ■ RT) 

70 C 

UO IF (SLOPEP .BQ. 0.0) 60 TO 120 

7INT «•' (RT - X) ! SLOPEP + Z 
XINT - RT 

IF (ZINT .6E. 0.0 .AND. ZINT .LE. THRLEN) GO TO 300 

75 C 

C TEST FOR INTERSECTIONS ALONG THE SPACE CRAFT SURFACE BETWEEN (0, RT) 
C AND (0# RBOUNO). (EQUATION OF BOUNDARY LINEI Z ■ 0.0) 

C 

120 IF (SLOPEP .LT, -l.OE+10) GO TO 130 

80 tint ■ 0.0 

XINT K X « 7 ♦ SLOPEP 

IF (XINT .GE. RT .AND. XINT .LE. RBOUNO) GO TO 300 
C 

C TEST FOR INTERSECTIONS ALONG RBOUNO BETWEEN THE SPACE CRAFT WALL 
85 C (0# RBOUNO) AND (7B0UND# RBOUNO). (EQUATION OF BOUNDARY LINEi 
C X - RBOUNO) 

C 

130 IF (SLOPEP .EO. 0.0) GO TO 400 
ZINT • (RBOUND - X) 7 SLOPEP + Z 
QO XINT ■ RBOUND 

IF (ZINT .GF. 0.0 .AND. ZINT .LE. Z80UN0) GO TO 300 
C 

C TEST FOR INTERSECTIONS ALONG THE BEAN EDGE BETWEEN THE THRUSTER 
C (THRLEN# RB) AND (ZBOUND# RB). (EQUATION OF BOUNDARY LINEI X - RB) 
95 C 

ZINT - (RB - X> / SLOPFP + Z 
XINT « RB 

IF (ZINT .GF. THRLEN .AND. ZINT .LE. ZBOUND) GO TO 300 
GO TO 400 

100 C 

C CALCULATIONS FOR RIGHT. 

C - 

c 

C TEST FOR INTERSECTIONS ALONG ZBOUND BETWEEN THE BEAM EDGE (ZBOUND# 
105 C RB) AND (ZBOUND# RBOUND). THE Z COMPONENT OF THE TOTAL VELOCITY 

e CHECKED TO DETERMINE THE DIRECTION OF PATH PROPAGATION, (EQUATION 
C OF BOUNDARY LINE» Z » ZBOUND) 

C 

200 IF (SLOPEP .LT. -l.OE+10) GO TO 210 
UO ZINT • ZBOUND 

XINT •• (ZBOUND - Z) ♦ SLOPEP + X 

IF (XINT .GE. PB .AND. XINT .LE. RBOUND) GO TO 300 
210 IF (VELTZ(D) 220# 400# 240 
C 
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C TEST FPP INTERSECTIONS ALONG RBOUND BETWEEN THE ^PACE CRAFT SURFACE 
C (0# RBOUNO) AND (ZBOUNO» RBOUND). (EQUATION OF BOUNDARY LINBI 
C X " RBOUND) 

C 

?;20 ZINT ■ (RBOUND - X) / SLOPEP + I 
XINT • RBOUND 

IF (ZINT .GE. 0.0 .AND. ZlNT »LE. ZSOUNO) GO TO 300 
GO TO 400 
C 

C TEST FOR INTERSECTIONS ALONG BEAN EDGE BETWEEN END OF THRUSTER 
C (THRLEN# RB) AND (ZROUNO# RB), (EQUATION OF BOUNDARY LINEI X ■ RB) 
C 

240 ZINT ■ (RB « X) / SLOPEP + Z 
XINT - RP 

IF (ZINT .GE. THRLFN .AND. ZINT .LE. ZBOUND) GO TO 300 
GO TO 400 
C 

C CALCULATE THE DISTANCE TO THE BOUNDARY. 

C 

300 DS • SORT ((Z - ZINT) ♦♦ 2 + (X - XINT) ♦♦ 2) 

IF (DS *GT. 0.25) GO TO 402 
RETURN 
C 

C ERROR CONDITIONS +♦♦*♦♦♦•> 

C 

C - ERROR 410 - NO "GOOD” INTERSECTION POINTS FOUND BETWEEN 

C BOUNDARIES AND CURRENT PATH USING DS (FATAL). 

C 

C - ERROR 412 - OS UNUSUALLY LARGE (NON-FATAL). 

C 

400 WRITE (inUT. 410) X, Z. SLOPEP, XINT, ZINT, I, VELTZ(I), LP 
60 TO 408 

402 IF (ICLERR .EQ. 1) GO TO 409 
IF (LR .EQ. 1) GO TO 409 

WRITE (IDUT, 412) LR, I, OS, Z, X, SLOPEP, XINT, ZINT 
GO TO 409 

408 OS » 1.0E+20 

409 RETURN 
C 

C ERROR CONDITION FORMATS. 

C 

410 FORMAT ( Z,11X,23H+»4+>«'^ ERROR 410 /, /, IIX, 

1 30HDS UNUSUALLY LARGE (NON-FATAL ),/, IIX, 

2 34HCALLE0 FROM FUNCTION SUBROUTINE DS,/,11X, 

3 3HX -,e9.3,5H 7 -,E9.3,10H SLOPEP ■,E9.3,8H XINT «,E9.3, 

4 8H ZINT -,F9.3,8H VELTZ(,I5,4H ) -,E9.3,6H LR -,I3) 

412 FORMAT (/,11X,23H****** ERROR 412 ♦*♦♦**,/,/, IIX, 

1 30HDS UNUSUALLY LARGE (NON-FATAL ),/, IIX, 

2 34HCALLE0 FROM FUNCTION SUBROUTINE DS,/,11X, 

3 4HLP -,I3,5H I -,I5,6H DS ",E9.3,5H Z -,E9.3,5H X -, 

4 E9.3,10H SLOPFP -,E9.3,8H XINT -,E9.3,BH ZINT -,E9.3) 

END 
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SUBROUTINE VRSPtT 
C 

C— VRSPIT USES VERSATEC PLOTTER TO PLOT ARPAVS XH AND YU 
C 

CHD THIS ROUTINE TS SITE DEPENDENT, IT PLOTS THE CONTENTS OF ARRAYS 
CWD XION AND ZION FROH POPE MEMORY, 

C 

C BLANK COMMON FOR LARGE ARRAYS 

COMMON nON( 4 l#l 5 U#XION( 4 U 151 ),VZU 5 X)»VXa 51 ) 
l»NIP( 4 l)» 0 N( 41 rl 5 lUDNH 42 UISTAT( 4 U 
COMMON / 10 t IN#I 0 UT#INF 0 { 14 )#KEY#ICLPLT>ICLWRT# mtL( 28 ) 

1 #IPATHS>IW,IF 1 ( 2 UIF 2 { 4 »#ICLERP 

C 0 HH 0 N/PARAM/N,NUHinN»NUMIT^RB#RfiDUND<RT#TBL 0 UT# 8 HCUR#UTU# 

1 TELIN,THRLEN>UMSION»VELBQH»ZBOUND»UjIR>PI»BK,Q,R 095 N 
1 # 0 N 06 #CEXSFC,TTHNEU»TIHE#TIHEMU»XVELNU# 2 VELMU»NSTAG# 

5 NSTGMU#NTOTST#PinV? 

DIMENSION SAV7(2)f SAVX(2)#D/a51),0XU51) 

DATA ZAXLN,XAXLN,INC,UNTYP#ISYM / 9 , 0 # 7 , 0 / 1 #+ 0 ^ 1 / 

DATA INTRZO/ 

CW USE INTR TO COUNT ENTRY NUMBER AND AVOID REINITIALIZING 
C— FIRST ENTRY-- SET UP THE SYSTEM^ SCALE, AXES AND TITLE IF DESIRED 
INTR-INTR +1 
IF(KEY.EQ, 0 ) GOTO 3 
IFdNTR.EQ,!) CALL PLOTS( 0 ,, 0 ,> 0 , ) 

C— SET ORIGIN OF PLOT 

CALL PLOT a.il./- 3 ) 

CALL SETHScm 
SAVZdl-O. 

SAVX ? 

SAVZ 12 )-ZR 0 UND/ZAXLN 

SAVX( 2 >-RB 0 UNP/XAXLN 

ZMAX-ZAXLN>I'SAV7(?' 

CW WORD-SIZE DEPENDENT AREA? IWl, IW 2 ARE FLAGS IN ITITL 
IW 1 -IW+ 1 
IW 2 -IW 1 + IW /2 

CALL AXIS( 0 ., 0 .,ITrTL(IWl),- 40 ,ZAXLN, 0 ,,SAV 2 (l)»SAVZ( 2 )) 

CALL AXIS( 0 .,XAY!.N, 1 H , 1 , ZAXLN, 0 .,SAVZ ( 1 ) , SAVZ ( 2 » ) 

CALL AXIS{ 0 ., 0 .,mTL(IW 2 )i 40 ,XAXLN# 90 .,SAVX(l»#SAVX( 2 U 
CALL SYMBaL( 0 .# 8 , 0 , 0 . 14 iITITLa»# 0 .» 80 ) 

DO 00 J- 2 #NUMI 0 N 
IFtZIQNU/D.GT.ZMAX) GO TO B 2 
80 CONTINUF 
82 NUHPT-iJ-l 

IF(K£Y, 6 T.O) NUHPT-NUHION^l 
DO 150 J« 1 #NUHPT 
NPTS-ISTAT(J) 

IF(ISTAT(J).E 0 , 0 ) NPTS-NIP(J) 

DO 100 NH- 1 ,NPTS 
OZ(NM)-ZION(J,NM) 

OX(NM>«XIDN<J»NH) 

100 CONTINUE 

DO 120 1 - 1,2 
NPT-NPTS+I 
DZJNPT)-SAVZm 
120 OX(NPT)-SAVX(n 

IF(J.E 0 .NUMI 0 N+ 1 ) lSYM -0 

CALL LINE CDZ,DX,NPTS, INC, LINTYP,ISYM) 
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150 CONTINUE 

C— DRAM SCHEMATIC OF THRUSTER AND BEAM, 

60 nxm»RT 

DX(2)«RT 

omi»o. 

DX(3)-0# 

DZ(2)-THRLEN 
65 DZ(3)-THRIEN 

0Z(AJ»THRLEN 
DX<A)*Rfl 
0X(5)-RB 
0Z(5)«ZB0UND 
70 NPTS«5 

00 200 I*l»2 
NPT-NPTSfl 
DZ(NPT)-SAVZm 
200 DX(NPT)«SAVX(n 
75 CALL UNE<DZ*OX#NPTSrl^O#0) 

CW FINISH THIS PLOT ANO 60 BACK FOP MORE 
CALL PL0T(0. ♦0.1-909) 

return 

CW-TERMINATE ALL PiaTriNG*i^>RELEASE OUTPUT TO VERSATEC PLOTTER 
BO 3 IF UNTR ,Lts RETURN 

CALL PL0T(0.»0,|t +999) 

INTR • 0 

RETURN 

END 
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1 SUrMlRUTIHE VRSf>l 

C 

C— VRS?l, USES VEBSATPC PLOTTER TO PLOT ARRAYS X<l ANO Y(» 

C 

5 CWO THiS ROUTINE IS SITE DEPENDENT. IT PLOTS THE CONTENTS OF THE 
CWO ARRAYS XION AND ^lOH WHICH HAVE REEN PEAD FROM THE EXTERNAL 
CWO FILE PATHS IN TERMS OF TRIPLETS. 

C 

C BLANK COMMON FOR LAPGF ARRAYS 
10 COMMON nON{AX,lBlbXION<Al#151J»Via51bVX<l51) 

l*NIPUl)#DNf Al|l51)iONHA2)#ISTAT(Al) 

COHHOH / 10 / IN>I0UT|lNFaa4)#KEYiICLPLT#iCLWRT#ITITL<28) 

I >IPATHS>IW#IFU2),IF2H>>ICLERP 

COMHDN/PARAM/N#NUHlON,NUHlT#RB#RBOUND»RT#TELDUT#BHCUR#UTUi 
15 I TBUN,THRLFN,UMSlONiVELflQH»ZBQUND»IL»IR»Pl#BK»Q#R8«?5H 

1 #ONOB»CEXSEC»TTHHEU»VTHE,TIMEMU>XVELMU#ZyELMU^NSTAG^ 

3 NSTGMU»NTOTST#PIOV? 

DIMENSION SAVZI2)>SAVXt2)>0Z(15U#DXU5U^0<lSll 
DATA ZAXLN#XAXLN,INC»LTNTYP,UYM /B.0»6.5f 1»430^1Z 
20 DATA INTRZO/ 

CW USE INTP TO COUNT ENTRY NUMBER AND AVOID PEIHITIALIZIHG 
C— FIRST ENTRY-^SFT UP THF SYSTEM# SCALE# AXES ANO TITLE IF DESIRED 
CWO U « 23 Bl» MUST REWIND FILE BEFORE READING 
REWIND IPATHS 
25 INTR"INTP+1 

IF(KEY.EQ.O) GOTO 3 
IFCINTR.EO.X) CALL PLOTSl0.#O.#O,) 

C— SET ORIGIN OF PLOT 

CALL PLOT (l.#l.#-3> 

30 SAVZ{ii**Pr. 

SAVXUHO. 

SAVZ(2)'»0.92/7AXLN 

SAVX(2)-R80UNn/XAXLN 

ZMAX-ZAXLN4SAV7(2) 

35 CW WORD-SIZE DEPENDENT APEAl IWl# IWZ ARE FLAGS IN ITITL 

iwi-ru+ 1 

IW2«IWl lW/2 

CALL AXIS(0..0,#ITITUIWl»#-A0#ZAXLN#0.#SAVZ(l)#SAV7(2l ) 

CALL AXISI0.#XAXLN#1H # 1# ZAXLN# 0. #SAVZ ( 1) #S AVI(2 ) ) 

AO CALL AXIS(0.>0.#ITITLtTU2)#AO#XAXLN#qO.#SAVXU)#SAVX{2)| 

CALL SYHB0L(1.»B*0,0.1A,ITITL<1)#0.#80> 

CW— 0C2981— HOD FOP READING# PLOTTING PATHS FILE 
WPITE{I0UT»2U NUMION#NUHIT 
NUHl - NUHtON 4-1 
A5 DO 55 JS- 1# NTOTST 

CWO THE FOLLOWING CODE COUNTS THE NUMBER OF TRAJECTORIES WRITTEN 
CWO TO FILE PATHS. 

ICOUNT • 0 

IF {JS -1) 5#15#S 

50 5 IFVAR - ((JS-1) ♦NUHIT +1) - ({JS-2)A10) 

DO 7 I»X#NUMIClN 
IF{NIP{I)-IFVAR) 7#7,9 
7 CONTINUE 

9 INITOO • T 

55 DO XX I « l#NUHION 

n • NUMl -I 

IF(NIP{II)-IFVAR) 1I»XX#X2 
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so 
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IT- 


90 
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no 


U CONTINUE 

12 LASTOO • n 

DO lA IJ •• INITOOi LASTOD 

IP lISTATUaj - SBABJ 14» 13> lA 

13 IPVAR - US NUHTT ♦ n *. H,I5 - 1) A lOJ 
IP (NlPtn - IPVAR) 6# 6# lA 

6 ICOUNT • ICOUNT + I 

14 CONTINUE 

LASTOO a lASTOO - INITOO « ICOUNT 
GO TO IS 

15 INITOO • 1 

lASTOO - NUMION 

IB CONTINUE 

DO 44 LASTOO 

REAOJIPATHS) I0N#NITER 
MRITFt 10UT»22) ION#NITER 

22 F0RMAT(26H +f+ VRSPl +4f ION NUHBER #I3n5»llH ITERATIONS/A) 
inteo-5 

21 F0RHAT(lH0n0i^»l3H +++ VRSPL +4+,I5#6H IONS #15|11H ITERATIONS ) 
DO 33 IT- IfNITER 

CWD 12/7/81 VX#V7 TO DX^DZ AND 0 TO OUT) TO HOLD DENSITIES 
REAOUPATHS) OX{lT)#DzaT)#D{lT) 

IFUOPUPATHS) ♦N0. 0,0) GO TO 70 
33 CONTINUE 
DO 37 I-l#2 
NRT - N1TER41 
0Z{MPT)-SAV7Cn 
57 ny(MPT)aSA\/y(T ) 

CALL LINE {OnDX,NITEP,lNC#UNTYP>INTEQ) 

44 CONTINUE 
55 CONTINUE 

C— ORAV SCHEMATIC OP THRUSTPR AND BEAM, 

70 DXm-RT 
DXi2)-PT 
OZ(1)-0. 

0X{3)-0, 

DZ(2)-THRLEN 

0Z(3)-THRLEN 

0Z(4)-THRLEN 

0X(4)-RB 

DXt5)-RS 

0Z{5)-ZB0UNP 

NPTS-5 

DO 200 I-l#2 
NPT-NPTS+I 
OZ{NPT)-SAVZm 
200 DX(NRT)-SAVX(n 

CALL LINE(DZ,nX,NO7SU»0.0) 

CW FINISH THIS PLOT AND GO PACK FOR HOPE 
CALL PLaTr0,»C,»-999) 

RETURN 

eW-TEPKiNATE ALL PLOtTING-»PELEASE OUTPUT TO VERSATEC PLOTTER 
3 IF UNTR ,LT. ?) RETURN 
CALL PLnT(0.>0.» +999) 

INTR - 0 

RETURN 

END 
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1 SUBROUTINE HOTW ( X,Y,r#LB#UNE»LAB) 

C 

C— SUB. PtOTW BY 0. R. WINDER# PHYSICS DEPT.# COLO. ST. UNIV. 

C PLOT ARRAYS Y VS X# PACH HAVIHC N POINTS# SELF SCALING. 

9 C 

C-- LB IS THE SYNPOL USFO FOR THE CURRENT GRAPH— SUGGEST HUHERICAL ORDER 
C— E.G.- FDR THE FIRST GRAPH# SET IN CALLING PROGRAHl LB-lHl 
C— THEN FOR 2ND ONE# RESET iTl LB«1H2# ETC. ONLY ONE CHARACTER# PLCASE. 
C«*4HE FIRST ENTRY IS CRITICAL— IT ESTABLISHES NUMBER OF POINTS# ALSO 
10 C— UPON FIRST ENTRY# MAX AND HIN VALUES ARE FOUND AND ARE USER LATER 
C— IF LATER CALLS INVOLVE POINTS OUTSIDE THESE LIMITS# THEY WILL BE 
C— TRUNCATED AND AN ERROR MESSAGE PRINTED ON THE PLOT FILE. 

C— TO SIGNAL I LAST GRAPH TO BE PLOTTED# PUT N-'<0. ALL GRAPHS WILL 
C— BE PLOTTED ON ONF SHEET VIA THE ARRAY LINE (103/60) . IF ONLY ONE 
15 C— GRAPH IS DESIRFO# SET N— N IN THE CALLING PROGRAM. 

C— NOTE ON SCALING* THIS EXPECTS PRINTER WITH 10 CHARS/INCH# 6 LINES/IN 
C— LA8(I6) CONTAINS TITLE (IN WORDS 1-8 )# X-AXIS LABEL (IN 9* 12) 
CP.-AND Y-*AXIS LABEL (IN 12-16). THE LAST CALL (N«0) DETERMINES LABELS 
C 

20 C— CALLING PROGRAM MUST SET UP L1NE(103#60) AND LABab) 

C 

C— NOTE THAT FORMATS ASSUME IWO-103# LN6»60. IF OTHERWISE# ADJUST THEM 
20 F0RMAT(8Hl PLOT# 22X#8A10#/#29X/8HX-AXISI #AA10#UH# Y-AXIS* # 

2 4A10 ) 

25 21 F0RMAT{20X#1HI»X0(10H,,..V....X)#2H.I) 

22 F0RMAT(//#7X#1AHPL0T— ON ENTRY#I3»12H WITH SYMBOL# IX# Al# /9X# 

2 52HTHE RANGE OF VALUES EXCEEDED THAT SET ON FIRST ENTRY#/#9X» 

3 31HCURRENT RANGES ARE* (ABSCISSA)#14X#10H(0RDINATE)#/»12X# 

4 13HFIRST ENTPY— #4X,4Eil .3#/#12X#12HTHrS ENTRY— #5X# 4EU.3# ///) 

30 23 F0RMAT(//7X#13HPL0T— ENTRY#I3#0H SYMBOL ,.Al»I9#14H POINTS# RANGE 

2 31HS OF ABSCISSA AND ORDINATE ARE »#/29X#4E11.3) 

24 F0RMAT(20X#103A1) 

25 F0RHAT(5X#1PB13.4»2X#103AI) 

26 F0RMAT(13X#U(IX#1P£9.1 )#///) 

35 C 

DIMENSION X(N)#Y(N)#LINE( 103/60) #LAB(16)#ZX( 11) 

DATA IBLNK#IB0RD#IW0#LNG#NTR/IH #1HI#103#98#0/ 

CW 7/27/81 REPLACED DATA ON I/O WITH NEXT LINE— RH PROBLEMS 
IOUTl-6 

40 C 

NTR-NTR + I 
IF(NTR.E0.1)NN*IABS(N) 

- MAX-MIN SECTION— FIRST ENTRY ONLY THE RANGE IS SET FOP X AND Y 
XSl-X(l) 

XLl-Xtl) 

YSl-Y(l) 

YL1»Y(1) 

DO 31 J-2/NN 

XS1“AMJN1(XS1#X(J)) 

XLl»AHAXi(XLl#X( J)) 

YSI"AHIN1(YS1»Y( J)i 
31 YLI-AMAX1(YL1»V(J)) 

IFINTR.GT.DGOTO 34 

xs-xsi 

XL*XL1 
YS«YS1 
YL-YLl 


45 


50 


55 

I 
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fiOTO 39 

C-—NOT 7HP FIRST PNTRY» Sn CHECK RANGES 
60 3A IFtXSl.LT.XS) GOTH 33 

IFtXLl.GT.XU GOTO 35 
IF(YSi.l,TtYS} GOTO 35 
IFIYU.GT.YU GOTO 35 
GOTO 39 

65 C--QUT OF RANGEw'^'WRITF MESSAGE AND SET A FLAG WITH NTR 
35 WRITEU0Un#23) NTR/LF#XS#XI.^YS#Yl.»XSl#XU?YSl»YU 
NTRp - NTR 
GOTO 55 

C—MRITE MESSAGE ASOUT THIS IGOOD) ENTRY 
70 39 yRlTEdOUTliRS) NTfi,U#N,XSl/XU#YSl>Ytl 

IF(NTR.NE.i) GOTO 55 

C--OONE WITH RANGING. NOW SCALING ON FIRST ENTRY ONLY 
ySCALE»(XL»-XS)/fFL 0 ATnW 0 - 2 n 
YSCALE»(YL-YSJ / /FL 0 AT(LNG- 1 H 
75 DO 70 4 » 1 HI 

70 IX(J)-XS ♦ FLOATt J-U^XSCALE*^ 10 . 

C— RLANK THE PLOT ARRAY t LINE ON THE FIRST ENTRY 
00 33 J» 1 #IVD 

00 33 K» 1 ,,LN 6 

BO 33 LINE(J>K).IBLNK 

C— flOROERS—LEFT AND RIGHT 
DO J» 1 #LNG 
LINEUf J)«IRnRD 
AA LINE(IWD;Jl«lBORn 
UK RR rnurrMitc 

^ <0 wirmr'r t «r-WM 

C-- FILL IN ARRAY LINE 
DO 57 I»l>NN 
IX-(Xm-XSWXSCALF+ 1.5 
IF(X(n,LT.XS)IX«»l 
90 IF(X(n.GT.XL)IX»IWD 

lY-Um-YSJ/YSCALF + .5 
IF(Y(I).LT.YS)IY<iO 

IF(Y(n.GT.YL)IY-LNG -1 
IY»LNG - lY 

95 57 LINEaX>IYJ-LR 

C— IF NOT THE LAST GRAPH» GO BACK FOR MORE 
IFdJ.GT.OGOTO 91 

C— MUST BE PLOTIMF — PRINT TITLE# AXES LABELS# AND TOP BORDER 
WRITEt I 0 UTI#R 0 )LAB 
lOO WRITE<I 0 UT 1 # 2 U 

C— PRINT Y-VALUES AND PLOT 
YVAL<»YL 4 YSCALE 
00 73 J«i#LNG #2 
YVAL*YVAL * YSCALE 

105 WRITE(I 0 UT 1 # 25 ) YVAL# <LINE(L#J)#L- 1 #IV 0 ) 

YVAL-YVAL - YSCALF 

73 WRITEIIQUTI.PA) (L INE ( L# J+1 ) # L-1# IWO) 

WRITE! I 0 UT 1 » 21 ) 

WRITE( 1 QUT 1 # 26 I ( 7 X(L) »L- 1 # 11 ) 

110 C— CLEAN UP THIS «ESS AND RETURN 
91 NTR"IABS(NTR) 

C— IF THIS IS THE LAST GRAPH# RE-INITIALIZE 
rF(N.LT.l)NTR -0 
RETURN 

115 end 
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1 SUaROUTINB LMolt 

CW HZtitn LNPLT HnOIPIgn for plasihi naw no fornai. parahfters 
C INPLT PREPARES ARRAYS X AND Y FOR PLOTTING VIA PLOTV 
C N IS NUHftFR DP FNTRIFS IN X AND Y ON FIRST ENTRY 
5 C N«0 SIGNALS LAST ENTRY, FOR ONLY 1 ENTRY# ENTER WITH -N 
OIHENSIDN LINFa.D''#60bLABa6>#Xa51?#T(l!JX) 

CW 7/2A/8X CHANGFO CDOF FOR PLASIM PLOTTING. «"-STP -INCREMENTS 
CW 7/2A/8X RLANK COMHON AND LABELED CDHHON FROM READER 
COMMON ZIOH?4X»iai)»XJON(4X#X5X)#VELn(l5X)#VELTX(X5X) 

10 X#NIP{4X)#DN(AXa5U>0NI(A2)#ISTAT(AX) 

COMMON/IO/ IN»inUT#INFaUA)#KEY»ICLPLT#ICLWRT»ITITL(28) 
a »IPATHS#IW#rFl{2)#IF2H)#lCLERR 

CDHMON/PARAM/ N#NUHiaN»NUMIT#RB#RBOUND#RT#TELOUT#BHCUR#UTILi 
3 TEL1N#THRLF.N#UMSI0N#VELB0H#ZB0UND#U#IR#PI#BK#Q#RB95N# 

15 ^ DNOB#CEySFr>TTHNEU#TrME#TIMEMU#XVELMU#ZVELMU#NSTAG» 

5 NST6HU#NT0TST#PI0V2 

DATA NTP#I0NSTPiITRSTP/0»2#2/ 

NTR-NTR+X 

IFINTR.RT.DGnTO 3 5 

20 CW 7/24/81 PUT LABELS IN LAB FROM ITITL 
IW2-2*rW 
00 33 J-l#IW2 
33 LAKJI-ITITLtJ) 

35 CONTINUE 

25 CW 7/24/8X SET UP DUMMy t$l FOR PLOTW USING XION»ZION 
CW 7/24/Pl IF OTHER PLOTS DESIRED, CHANGE NEXT LINES 
NIT-NUMIT+1 

CW 7/28/01 AD HOC SFT UP MAXHIN FOR PLOTW 
DO 39 ION-1, NMMinN 

30 NrTfi-NrP(TnNJ-.(NSTAG-l)>l‘(NUHIT-10> 

X(l) - XION (IDN,1) 

Z(l) - ZION noN,i) 

x(NiT) • xa» 

ZtNIT> - Z(l) 

35 DO 39 ITR-?,NITR 

XU) - AMrNl(XU),XION(ION,ITP)) 

ZU) - AHrNl(ZU),ZtONUON,ITR)) 
X(NIT)-AMAyi{X(NIT),XIONUON,ITR) ) 
Z(NIT)-AMAyi(nNIT),ZIQNUON,ITR) ) 

40 39 CONTINUE 

LB-IH 

CALL PLOTW(?»y,NIT,LB,LINE#LAB) 

DO 6A IDN-2,NUMI0N#T0NSTP 
DO 44 J-I,NIT 
45 X(J)>0,0 

44 Z(J)-0.0 

NITR-NIP(ION) 

LB«SHIFT{ION,+54) 

DO 55 ITR-l,NITR,TTRSTP 
50 X(ITR)-XION(IDN,ITR) 

55 ZaTR)-ZlQN{ION,ITR) 

CW 7/24/01 IF LAST ENTRY TO PLOTW, TELL IT SO WITH N-0 
IFnON,GE.NUMinN)NIT-0 
66 CALL PLOTW(Z»y,NIT,LB,LINE#LAP) 

§5 CW 7/24/Pl FINISHED WITH ONE PLOT# RESET FOR NEXT ONE 
NTR-0 
RETURN 
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